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Abstract 

A detailed numerical study is presented of the slow diffusion (Arnold diffusion) taking place 
around resonance crossings in nearly integrable Hamiltonian systems of three degrees of free- 
dom in the so-called 'Nekhoroshev regime'. The aim is to construct estimates regarding the 
speed of diffusion based on the numerical values of a truncated form of the so-called remainder 
of a normalized Hamiltonian function, and to compare them with the outcomes of direct numer- 
ical experiments using ensembles of orbits. In this comparison we examine, one by one, the 
main steps of the so-called analytic and geometric parts of the Nekhoroshev theorem. Thus: i) 
we review and implement an algorithm |20] for Hamiltonian normalization in multiply resonant 
domains which is implemented as a computer program making calculations up to a high normal- 
ization order ii) We compute the dependence of the optimal normalization order on the small 
parameter e in a specific model and compare the result with theoretical estimates on this depen- 
dence, iii) We examine in detail the consequences of assuming simple convexity conditions for 
the unperturbed Hamiltonian on the geometry of the resonances and on the phase space structure 
around resonance crossings, iv) We discuss the dynamical mechanisms by which the remainder 
of the optimal Hamiltonian normal form drives the diffusion process. Through these steps, we 
are led to two main results: i) We construct in our concrete example a convenient set of variables, 
proposed first by Benettin and Gallavotti [4], in which the phenomenon of Arnold diffusion in 
doubly resonant domains can be clearly visualized, ii) We determine, by numerical fitting of our 
data the dependence of the local diffusion coefficient D on the size \\Ropt\\ of the optimal remain- 
der function, and we compare this with a heuristic argument based on the assumption of normal 
diffusion. We find a power law D cc \\Ropt\\^^^^''\ where the constant b has a small positive value 
depending also on the multiplicity of the resonance considered. 
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1. Introduction 

The study of diffusion in nearly-integrable Hamiltonian dynamical systems of the form 

H(I,<p) = HQiI) + £Hi(I,(p) (1) 

where (/, <p) are n-dimensional action - angle variables and e is a small parameter, constitutes a 
central problem in Hamiltonian dynamical systems theory, in view, in particular, of its multiple 
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applications in physics and astronomy (see for an introduction, the basic review paper 

ImI^ or 1.47.1 bTHi llsgll for recent advanced reviews emphasizing various aspects of this subject). It 
is a well established result that, if « > 2, and H satisfies appropriate convexity and analyticity 
conditions (see section 2 below), two distinct regimes characterize the laws of diffusion as a 
function of e: for e < e*, where e« is a threshold value, the onset of the so-called 'Nekhoroshev 
regime ' takes place EHnHlHl Qilllii. In this case, the Nekhoroshev theorem provides an 
0[exp(-(e,/e)'^)] upper bound for the speed of diffusion. The exponent c depends on the number 
of degrees of freedom «, while its precise value in local domains of the action space depends also 
on the multiplicity of the resonance conditions holding in such domains (see e.g. i4ll ll54llll48ll ). 
Furthermore, the mechanism of diffusion caused by transition chains, as demonstrated in one 
special example by Arnold (see also Ii561 ). is conjectured to hold in more general systems of 
the form ([ij (e.g. ll46ll : note, however, that no formal proof of this fact has been given to date). 
On the other hand, for e > e», the diffusion is driven mainly by the mechanism of resonance 
overlap ll55|]|fl3|] In this case, one expects a power-law dependence of the speed of diffusion 
on e (see e.g. (|9|]; a power law is also found in the case of the so-called 'Fast Arnold diffusion' 
iP). 

The diffusion in weakly chaotic systems has been a subject also of extensive numerical stud- 
ies over the last three decades (some indicative references are ll36llll59ll l37ll 151 llssil 16lll25llll7ll'). 
A detailed study, however, of the very slow diffusion characterizing the 'Nekhoroshev regime' 
has become possible only in recent years. In this respect, we note in particular the series of in- 
structive works II22II II38IIII23IIII33IIII35II . where, using the so-called Fast Lyapunov Indicator (FLI; 
see II22II '). a method was found to depict the resonant structure of the action space in models of 
three degrees of freedom, or 4D and 6D symplectic mappings being in the Nekhoroshev regime 
1 22, 23, 35]. In llssi] . the mean-square spread in action space < A/^ > was measured as a function 
of the time f for orbits along the chaotic border of a simply-resonant domain (see section 2 for a 
precise definition). It was found that i) the local character of diffusion is normal, i.e. < >oc t, 
and ii) the diffusion coefficient D -< AJ^ > /t decreases with e faster than a power law. The 
exponential fit D oc exp(-(e(./e)"-^^ was given in a subsequent study (l^. The estimate obtained 
in 1,35.] . through interpolation over five orders of magnitude of the perturbation parameter, yields 
with with certainty the first digit of the exponent 0.2 . . ., but the errors in the interpolation make 
uncertain the second digit in both the above estimates. In [41], D was measured as a function 
of the separatrix splitting S of the asymptotic manifolds of simply unstable two-dimensional tori 
lying at the borders of simple resonances (see also |49|(4^). The measurement of S itself was 
based on employing the FLI. It was found that D cc S'', with p - 2.1 and p - 2.56 in two reso- 
nances of increasing order respectively. Finally, the laws of diffusion in systems violating one or 
more necessary conditions of the Nekhoroshev theorem were investigated in 1 34]|42l, leading to 
a number of interesting results regarding the dynamical consequences of such violations. 

The motivation for the present study stems primarily from the results reported in refs L38l ll23ll 
llisll 1 41 ], and it can be described as follows. The results obtained so far are very satisfactory from 
the numerical point of view. They require, however, computations involving large ensembles of 
orbits and integration times of the order of billions, or even trillions of periods. On the other 
hand, we can remark that, in principle, the analytical methods involved in the main theories of 
chaotic diffusion lend themselves also conveniently to getting quantitative predictions regarding 
the value of the diffusion coefficient, or the scaling laws of diffusion, in general, in the weakly 
chaotic regime. For such a goal, however, to be accomplished, it is required that one should be 
able to carry on expansions of certain quantities up to a very high order in the small parameter 
e (usually with the aid of a computer). This fact is explicit in Nekhoroshev theory, where one 
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needs to reach an expansion order high enough for the asymptotic behavior of the perturbation 
series to show up. This has been realized in studies seeking to determine the range (in the small 
parameter value) and/or the conditions of applicability of Nekhoroshev theory, or, finally, the 
domain of practical stability for motions in simple physical systems or models inspired mainly 
from Solar System dynamics (see e.g. l2M\Mi L48J LISJ LiSJ LIS] L43J152JL31J). These studies 
notwithstanding, the question of central interest in the present paper, namely how to obtain rel- 
evant quantitative estimates of the local value of the diffusion coefficient D in resonant domains 
(of various multiplicities) of the action space via high order expansions of perturbation theory, 
remains, to our knowledge, largely unexplored. 

Regarding now this last question, it should be noted that the formal analytical apparatus of 
Nekhoroshev theory, entailing the construction of a normal form in local domains covering the 
action space of systems of the form ([T]i, aims to transform the original Hamiltonian into one 
in new canonical variables resuming the form Htrcmsformed - Z + R, where Z, the normal form, 
corresponds to a simple dynamics, while R, the remainder, induces a perturbation to this dynam- 
ics. The so-called 'geometric part' of Nekhoroshev the theorem ensures that, despite allowing in 
general for chaotic motions, the flow under a multiply-resonant normal form alone would imply 
perpetual confinement of all chaotic orbits in balls of radius (9(e'^^) in the action space. Nev- 
ertheless, this picture is altered due to the effects of the remainder which eventually causes the 
orbits to diffuse away of their initial (9(e''^) domain. Now, via a sequence of hamiltonian nor- 
malization steps we find that there is an optimal order at which the size of the remainder becomes 
exponentially small in a power of 1 /e. This, in turn, implies an exponentially small semi-analytic 
upper bound of the value of the diffusion coefficient D. Unfortunately, such a bound turns usually 
to be very unrealistic, as it overestimates by a large factor the true value of D (or, equivalently, it 
underestimates the time of practical stability). We are thus led to conclude that, whereas the re- 
mainder R constitutes a quantity of primary interest in quantitative applications of Nekhoroshev 
theory, the precise relation between R and D is apparently very different from what upper bound 
estimates would suggest. Instead, a detailed analysis of the effects of the remainder on dynamics 
appears to be necessary in order to formulate a more precise theory of the relation between R and 
D. 

In the sequel, we present such an analysis in systems of three degrees of freedom. In this 
analysis, we still have to rely on an assumption for which numerical indications are available, 
namely that the local character of diffusion in sufficiently small domains of the action space is 
'normal', that is, the mean square spread of the actions of the chaotic orbits grows linearly with 
time (there are indications that global diffusion, which concerns ensembles or orbits diffusing 
in a substantial part of the Arnold web over much longer timescales, could also be described as 
'normal' (see ll33ll ): however, the issue of the laws of global diffusion can only be hoped to tackle 
after the laws of local diffusion have been adequately understood). In the rest of our analysis, we 
proceed by expressing all quantities of interest in terms of the remainder function, which, in turn, 
is calculated in concrete examples by a well-defined algebraic procedure. Finally, we estimate 
via this analysis how D depends on the size ||/?o;«ll of the remainder at the optimal normalization 
order It should be noted that the idea that the stability properties of the orbits in nearly-integrable 
systems depend on the size of the optimal remainder is not new, but it is one permeating nearly 
all forms of canonical perturbation theory. The novel feature here, instead, is to use ||^„y,f|| not 
as an upper bound for D, but as a way to estimate D via examining the relation between the 
two quantities as determined by independent numerical experiments. One main prediction is that 
this relation is altered according to the multiplicity of resonance conditions holding in the action 
domain of interest. More concretely, we predict that the diffusion coefficient D scales with ||^op/|| 
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as a power-law D oc where /? ^ 2 in doubly-resonant domains, while p - 2(1 -i- b) in 

simply resonant domains, for some constant b > 0. A combination of theoretical arguments 
found in (1^1 n I, together with quantitative estimates on the relation between the size of the so- 
called separatrix splitting (see subsection 2.3.2) and the normal form remainder given in ll49ll . 
suggest b ^ 0.5, i.e. /? ^ 3 in simply resonant domains. This agrees with the numerical results 
obtained in a previous study lEoll . 

In [id], a computer-algebraic program was written in order to calculate the optimal normal 
form as well as the remainder function Rgp, at the optimal mormalization order in a case of 
simple resonance, employing the same Hamiltonian model as in [38]. This operation involved 
computing about 5x10^ Fourier coefficients, at a truncation order in Fourier space as high as 
K - 44. Comparing the computed size of \\Ropt\\ versus available numerical data on D from [38ll . 
the scaling D oc ||/?oy«|p was found by numerical fitting. In the present paper, after presenting 
some theoretical results, we make a similar numerical calculation as in [20] but in the case of a 
double resonance. In order to reach the optimal normalization, we had to extend all normal form 
calculations up to the Fourier order = 50 (8 x 10^ coefficients). We thus determined the size 
of the optimal remainder \\Ropt\\ for many different values of the small parameter e. In the same 
time, we computed the diffusion coefficient D for the same values of e by a purely numerical 
procedure involving runs of ensembles of chaotic orbits (see section 3). Finally, we made two 
independent numerical comparisons of the relation between D and ||^oy„||. The latter yield the 
power laws D oc ' and D oc ||7?op/||^'' respectively. This essentially conffims that p is 

close to 2 in doubly resonant domains, albeit with a small noticeable difference even in this case, 
which probably requires a more precise theory to interpret. 

Besides the above computation, our analysis using high order normal forms resulted in a 
relevant result regarding the possibility to visuahze how the phenomenon of Arnold diffusion 
proceeds locally, within a doubly-resonant domain, by materializing the computation of a con- 
venient set of variables helping to this purpose, that were proposed in the work [4J. We note that 
numerical evidence for Arnold diffusion of orbits entering from simple to double resonances was 
presented in ^3^. Here, we provide a detailed topological description of this phenomenon. The 
whole computation consists of: i) computing a set of resonant canonical action-angle variables 
via a sequence of Lie canonical transformations, ii) taking a 2D Poincare surface of section of the 
doubly-resonant normal form dynamics (which represents a system of two degrees of freedom), 
and (more importantly) iii) using the energy Ez of the normal form as the third variable, showing 
the effect of Arnold diffusion. According to theory, the value of Ez changes exponentially slowly 
in time due to the effect of the remainder In the sequel we refer to this phenomenon as 'drift', 
although in reality it means that a number of quantities can be characterized as undergoing ran- 
dom walk during the whole diffusion process. Besides setting the timescale of diffusion, the 
drift can be viewed also as the source of a dynamical phenomenon, namely the communication 
between chaotic domains that would be otherwise isolated under the doubly-resonant normal 
form hamiltonian flow. We show in a true example the excursion of a chaotic orbit within the 
doubly-resonant domain as it appears in the above proposed set of variables. We thus identify 
a sequence of chaotic transitions of such an orbit from one resonant domain to another In fact, 
in each transition the orbit bypasses the barriers imposed by normal form dynamics via a 'third 
dimension', i.e. the slowly drifting value of Ez- We finally argue that, besides their practical 
utility, such illustrations are also suggestive of the geometric structure underlying the asymptotic 
manifolds of lower-dimensional tori filling the phase space in the domain of a double resonance. 
These manifolds are important, because, following the spirit of Arnold's original work [2|, it 
has been widely conjectured that their heteroclinic intersections constitute a primary cause of 
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Arnold diffusion. Of course, proving this fact represents a well known important open problem 
of dynamical systems' theory. 

The structure of the paper is as follows: Section 2 presents the theory, focusing on the normal 
form algorithm, multiply-resonant dynamics, effect of the remainder, and, finally, on the relation 
between D and ||/?opf||. We describe in some length all necessary theoretical steps in order to 
render the paper as self-contained as possible. Section 3 then passes to the numerical results. We 
present i) the results from the normal form computer-algebraic construction, ii) the visualization 
of Arnold diffusion using appropriate variables based on the normal form computation, iii) the 
numerical calculation of the diffusion coefficient D, and, finally iv) the comparison of D with 
||7?op,||. Section 4 summarizes the main conclusions of the present study. 

2. Theory 

Most statements made in subsections 2.1 and 2.2 below, regarding the properties of the 
Hamiltonian models considered as well as the algorithm by which we perform Hamiltonian nor- 
malization, are applicable to systems of an arbitrary number of degrees of freedom. In order, 
however, to be consistent with the rest of the paper, we use everywhere a notation referring to 
systems of three degrees of freedom. On the other hand, the analysis of subsection 2.3 applies 
to the study of the diffusion in doubly or simply resonant domains. In systems of three degrees 
of freedom, the latter represent the only possible multiplicities of a resonance condition, while 
in systems of more than three degrees of freedom there are also cases of intermediate resonance 
multiplicities between one and the maximal. The latter's study, nevertheless, is well beyond our 
present computational capacity, and thus it is left as an open problem. 

2.1. Definitions 

We consider three degrees of freedom systems of the form (H), where H satisfies the following 
analyticity and convexity conditions: 

i) Analyticity: H is assumed to be an analytic function in a complexified domain of its argu- 
ments. Namely, we assume that there is an open domain J c R-' and a positive number p such 
that for all points /» = {Iu,h*,h*) ^ I and all complex quantities I[ = Ii - lit satisfying the 
inequalities \r. \ < p, the function Hq admits a convergent Taylor expansion 

3 3 ^ 

Ho ^ Ho. + CO,- 1' + J] J] 2^iJ*I'I'j + ■■■ (2) 

1=1 M 

where = VjHoih), and M,y» are the entries of the Hessian matrix of Hq at Furthermore, 
we assume that there is a positive constant cr such that for all I € I, Hi admits an absolutely 
convergent Fourier expansion 

Hi^Yj^ik{I)txm-<t') (3) 

k 

in a domain where all three angles satisfy < Re{(pi) < 2n, |/m(0,)| < cr. By the Fourier theorem 
(see e.g. |30]), this condition implies that the coefficients hkiF) decay exponentially with the 
L'-modulus \k\ = -H \k2\ + 1^31, that is, there is a positive constant A such that the bound 

\h{r)\<Atx^{-\k\(T) (4) 
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holds for all k e We finally assume that all coefficients admit Taylor expansions with 
respect to /, 

1 ^ ^ 

(where htjjt. are the entries of the Hessian matrix of hk{I) at /,), which are convergent in the same 
union of domains as for Hq. 

ii) Convexity: For the Hessian matrix M,, which is real symmetric, we assume a simple quasi- 
convexity condition, namely that for all h e I either two of the (real) eigenvalues of M, have 
the same sign and one is equal to zero, or all three eigenvalues have the same sign. Furthermore, 
we define two constants: 

ju„„„ = min{|jUj|), yU„„.,- = max{|jt/j|) (6) 

where j is a label of only non-zero eigenvalues jij of M,, i.e. j - 1,2 or j - 1, 2, 3 if there are 
two or three non-zero eigenvalues respectively. 

As will be discussed in detail in subsection 2.3, the quasi-convexity condition is essential, 
since it introduces a confinement of the orbits for exponentially long times on a surface arising 
from the condition of preservation of the energy (see [4]). 

We now give some definitions allowing to characterize resonant dynamics. 

A resonant manifold Kk associated with a non-zero wavevector k with co-prime integer com- 
ponents k = (k\,k2, k^) is the two-dimensional locus defined by 

= {/€ J:^lW,(/) + fc2W2(/) + ^3W3(/) = 0) , (7) 

where w/(/) - dHo/dli. 

Let It G J be such that all three frequencies w,(/t), / - 1, 2, 3 are different from zero. We 
now distinguish the following three cases: 

i) Non- resonance: no resonant manifold !Rt contains 

ii) Simple resonance: one resonant manifold contains 

iii) Double resonance: more than one resonant manifolds contain In the latter case, it 
is possible to choose two linearly independent vectors k^^\k^^^ such that all resonant manifolds 
"Rk containing /, are labeled by vectors k which are linear combinations of the chosen vectors 
k^^\k'^^ with rational coefficients. The intersection of these manifolds forms a one-dimensional 
resonant junction. A doubly-resonant point /, always corresponds to the intersection of a reso- 
nant junction with a constant energy surface Ho(It) - E. 

In the above definitions, resonant manifolds of all possible wavevectors k have been con- 
sidered. It is well known, however, that in normal form theory a natural truncation limit \k\ < K 
arises in Fourier space (see below). Accounting for this possibility, we call a point /, e J i) 
non-resonant, ii) simply resonant, or iii) doubly resonant with respect to a K-truncation, if the 
number of resonant manifolds "R^ with \k\ < K passing through /» are i) zero, ii) one and iii) more 
than one respectively. 

Finally, it will be convenient to introduce a definition concerning open domains in I. Let 
'Wi^^B be a ball of radius B around one point /» in I. If Hq satisfies convexity conditions as 
assumed above, for B small whatsoever the domain ^i,,b is crossed by a dense set of resonant 
manifolds Hk- However, for any fixed value of the positive integer K, only a finite subset of 
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the manifolds Hk satisfy \k\ < K. The domain ^i,^b is then called: i) non-resonant, ii) simply- 
resonant, and iii) doubly-resonant with respect to the ^T-truncation if /» is, respectively, non- 
resonant, simply-resonant or doubly-resonant, and no other resonant manifolds with \k\ < K 
cross 'Wj^^B except for the ones passing through 

2.2. Normal form construction 

All our estimates on the speed of diffusion are based on an appropriate normal form con- 
struction. In this, we adopt the method exposed in detail in [20], which lends itself conveniently 
to i) developing a computer-algebraic program, and ii) deriving analytical estimates on the size 
of various quantities appearing in the course of Hamiltonian normahzation. The main elements 
of this method are: 

Expansion centers. The action space can be covered by domains 'W/. b, centered around 
points /, which serve as expansion centers of both the original Hamiltonian and the normal 
form. We choose the points to belong to the set of all doubly-resonant points of I, denoted 
by £), and by setting B as of order (9(e'^^). The covering is possible because D is dense in 
I. A normal form construction as done below is valid within one domain "Wj^^b (this is essen- 
tially the same starting point as in Lochak's |45] analytic construction leading to a proof of the 
Nekhoroshev theorem). A crucial remark is that the characterization of dynamics within 'Wj^^b 
as non resonant, simply resonant, or doubly resonant depends on e. This is because, as shown 
below, the optimal normal form truncation order K = Kgpi in Fourier space depends on the value 
of e. Furthermore, for a given value of K, the set D can be decomposed in three disjoint sets 
D = £)o,a'U£)i,a:U£)2,a:, containing all non-resonant, simply resonant and doubly resonant points 
respectively with respect to the /T-truncation. Thus, the characterization of resonant dynamics 
within depends on whether, according to the value of K, /, belongs to £)o_a:, Di^k, or D2.k- 

Resonant module: Let /. be a point of D and k'-^^ = (k\^\ k[^\ 4^'), /t^) = (kf\ kf\ kf) two 
linearly independent vectors such that fc*'' ■ coih) - for i - 1,2. More than one choices of A:*'' 
and A;*^-* are possible. In the sequel we choose A:*'' and k^^^ so that + \k^^\ is minimal. The 
vector m = (m\,m2, mj,) defined by 

k^k^-k!^k!^ (8) 

is parallel to the vector CL){h) since k ■ m - for all k satisfying k ■ cl){I,) - 0. If mi, m2, m3 are 
not co-prime integers, we re-define m by dividing the m, by their greatest common divisor The 
set 

M = [keZ^ -.k-m^o] (9) 

is hereafter called the resonant module associated with the point /, e £). The resonant module 
includes wavevectors k whose respective trigonometric terms exp(/A; ■ (p) are to be retained in the 
normal form. 

Action re-scaling: From now on we focus on the construction of the normal form in one 
specific domain W/,,b. It has been mentioned already that it is convenient to choose B as a 
quantity scaling proportionally to e'^^. The simplest way to accommodate such a choice is by 
introducing the following re-scaling of all action variables within 'W/,,^: 



Ji^e-''\li-Ii,)^e-"^i;, ( = 1,2,3 
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This re-scaling greatly simplifies the normal form algorithm, because it formally removes all 
terms besides linear in the actions from the kernel of the so-called homological equation (see 
below, or ll20ll for details) by which the normalizing generating functions are determined. Eq.lfTOt 
does not define a canonical transformation. However, the correct equations of motion in the 



variables {J, (p) are produced by the Hamiltonian function h{J, 
(neglecting a constant) 



^'^H(h + €^'^J,<f>), i.e. 



3 3 



h(J,ct>) = co.-J + e''^Y,Tu\^U*JiJj + --- 



(11) 



'=1 ;=i 



+ e''^ ^'k' + f '^'V/./i, ■ -/ + 5 X Z ^kJjJiJj + ■■■ 



3 3 



;=1 j=l 



exp(/A; ■ 



where the first line in the above equation comes from the integrable part Hq of the original Hamil- 
tonian (Eq.(l2]i), while the second line comes from the perturbation H[ (Eq.Q) given the series 
expansion of the Fourier coefficients as in Eq.©. 

Book-keeping: We now split the Hamiltonian (ITTT i in parts of different order of smallness, 
which are to be normalized step by step. The function (fTTT i contains terms of various orders in 
the small parameter e'^^. However, the presence of a second 'small parameter' e °" is implied in 
( fTTT i by the exponential decay of all Fourier coefficients hkt, hkjjt, etc., due to Eq.(|4ll (see jsoll 
pp. 90-91 for a thorough exposition of the role of this small parameter in Nekhoroshev theory). 



We take both parameters into account by introducing an integer K' such that e 
setting: 

1 



la- 



Me) 



e'^^, i.e. by 



(12) 



Using K' , the Hamiltonian ( fTTT i can be split in groups of practically the same order of smallness. 
This is realized by artificially introducing a 'book-keeping' coefficient AP in front of each term 
in (fTTI) . whose numerical value is set equal to unity at the end of the calculation. Furthermore, 
for a term of the form e^^'^fi J) exp(/A: ■ 0) we set p - {\k\IK'~\ + p. 

Regarding the above 'book-keeping' process it is worth noting the following: i) This way of 
splitting the Hamiltonian in different orders of smallness results in a finite number of terms ap- 
pearing in every power of A. ii) This technique is suggested already by Poincare 115 311 and Arnold 
Uj]. In fact, the dependence of K' on e is weak, since it is logarithmic, so that an alternative 
choice to the 'ansatz' (fTSl i is to set K' - const ~ l/cr. In fact, according to Giorgilli [30] this is 
an optimal choice, iii) Since, at every normalization order, we have a reduction of the analyticity 
domain, one could consider re-defining K' at every normalization step. However, this is hardly 
tractable from an algorithmic point of view. Instead, keeping K' constant at all normalization 
orders should be viewed as a rule indicating the sequence by which the various terms in the 
Hamiltonian are normalized, i.e., the terms or order A'' are normalized in the r-th step. Albeit not 
necessarily optimal regarding the grouping of the terms according to their size, this rule proves 
simple to implement and sufficient in practice. 

Returning to the form of the Hamiltonian after introducing the book-keeping factor A, the 
Hamiltonian reads: 



h 



1 



(13) 



J/2 3 3 \ 

+ ^2™\v,./.,.7+i^^«/'^'ii-y y ^^^^^^^^^^ . 



Setting Zo = w« ■ 7, the Hamiltonian ( fTll ) resumes the form 

oo 
.5=1 

where i) the superscript (0) denotes zeroth-step of the normaHzation procedure (- original Hamil- 
tonian), ii) the exponent of A in different terms keeps track of their true order of smallness, and 
iii) the functions Hf^ are of the form 

^(0) ^f^e^n" '"g" ' ^0)(y) exp(/fe . 4>) (15) 

fi=l \k\=K'(s-,,) 

where H^^^iJ) are polynomials containing terms of degree yu - 1 or // in the action variables J. 
Precisely, we have: 



if > 0, or 



i/W(7)=y y y — - — j^ji'jt 



if;t = 0. 



Hamiltonian normalization: We use the algorithm of composition of Lie series in order 
to perform the Hamiltonian normalization. Let us recall that the purpose of the normaliza- 
tion is to introduce a sequence of canonical transformations {J,<p) = (7*"', 0**^^) (7''\(^*'') 
— > {J^^\(p^-^) — > ... so that the Hamiltonian expressed as a function of the new variables al- 
lows one to more easily identify the main features of dynamics. After r normalization steps, the 
old variables {J,<p) = (y'*'\0*"^) are expressed in terms of the new variables (j'-''\<p^''''), and the 
Hamiltonian H^''\j^''\ <p^'''>) = h(J(f '\ <^(''>), 4>( J^''\ <*<'>)) takes the form 

H^''\ f \ </.*'■>) = Z^^H j^'K .^C); a, e) + R^''\/'\ A, e) . (16) 

The terms Z^''\j^''\ (p^''^; A, e) and R^''\j^'''\ <p^''^; A, e) are called the normal form and the remainder 
respectively. The normal form is a finite expression which contains terms up to order r in the 
book-keeping constant A, while the remainder is a series containing terms of order A''*^ and 
beyond. The mathematical structure of the normal form term Z'''' is such as to imply an easily 
identifiable dynamics in the variables ( J^''\ 0*'') (e.g. an oscillator or pendulum dynamics). On 
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the other hand, the remainder is a convergent series in a restriction of the domain of analyticity of 
the original Hamiltonian, which represents a perturbation with respect to the Hamiltonian flow of 
Z^''\ An optimal normalization order ropt exists (see below) where the process must be stopped. 
The Hamiltonian normalization is implemented step-by-step by the recursive equation: 

= exp(L^,.)i/*'-" (17) 

where is the r-th step Lie generating function and L^, = {■,Xr} is the Poisson bracket operator 
Both Z/*'^* and Xt- are functions of the variables The generating function is defined by 

the solution of the homological equation 

{cj, ■ f\xA + H'-r'\J^''\ 0^''') = (18) 

where H^r'^\j^''\ (jf''^) denotes all terms of //*'^'' which i) have a book-keeping coefficient /f in 
front, and ii) belong to the range of the operator {w, ■ J^''\ ■]. Given the definition of the resonant 
module M in Eq.©, one has the relation 

H^;-^^ ^ H^;-^^ - Zr (19) 

where H^^'"^^ are all the terms of //^'^"'* having a factor A'' , and Z, are the normal form terms of 
H^r ^^\ that is all the trigonometric terms whose wavevectors k belong to M. It follows immedi- 
ately that H^''^ has the form 

= Zo + Zi + ... + Z, + R^''^ (20) 

where all terms in the functions Z, have a factor A', while 7?^''' is a series in powers of A starting 
with terms of order /l'"^' . 



Optimal truncation: In the analytical part of the Nekhoroshev theory it is demonstrated that 
the whole normalization process has an asymptotic character Namely, i) the domain of conver- 
gence of the remainder series R^''^ shrinks as the normalization order r increases, and ii) the size 
of R^''\ where || ■ || is a properly defined norm in the space of trigonometric polynomials 
(see below), initially decreases, as r increases, up to an optimal order ropt beyond which ||/?*'^'|| 
increases with r. In the Nekhoroshev regime, one has ||Z*''"'"^|| >> Thus, stopping at ropt 

best unravels the dynamics, which is given essentially by the Hamiltonian flow of Z*''"'"' slightly 
perturbed by R^'''>p'\ The long term consequences of this perturbation, which determine the speed 
of diffusion, will be analyzed in subsection 2.3. 

The normal form Z*'^ = Zo + Zi + ...Zr contains trigonometric terms exp(/A: ■ 0) of order not 
greater than K - K'r - 1. Let ropt be the optimal normalization order It is well known that the 
dependence of ropt on e is given by an inverse power-law, namely 

ropt ~ . (21) 
The exponents 1/6, 1/4 and 1/2, referring to the non-resonant, simply resonant, and doubly 



resonant normal form constructions respectively, are found in 115411 . We emphasize that, while, 
due to the introduction of the book-keeping process, the algorithm of Hamiltonian normalization 
analyzed above is not technically identical with the usual normalization procedure used in the 
proof of the Nekhoroshev theorem (e.g. as in [54]), in practice we recover the estimate dJTJ, and 
the resulting exponents, both in the simply resonant case (see |'20'l) and in the doubly resonant 
case, as confirmed by numerical experiments in section 3 below. In particular, we find that since 
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the leading terms in the remainder are 0{A''"'"^^), the size of the remainder is of order Oie^''"'"'^^^^^), 
implying (viz.Eq.fO): 



i.e. the remainder is exponentially small in 1/e in accordance with the Nekhoroshev theorem. 
The Fourier order 



is hereafter called the optimal K-truncation order All the normal form terms of H^''"i"^ have 
Fourier orders satisfying \k\ < K„p,(e). 

2.3. Resonant normal form dynamics and the rate of diffusion 

We are now in a position to discuss the essence of all the previous definitions. The key point 
is to observe that, depending on the value of e, the same expansion point /, e 2) of the normal 
form construction turns to be either non-resonant, or simply or doubly resonant with respect 
to the optimal K-truncation. In particular, let fc^'^ and fc*^' be two linearly independent vectors 
of M such that for k e M one has \k\ > \k^^^\ > |^*''|. We then distinguish the following 
three regimes: i) < Kop,(e). Then, the point /, is doubly-resonant with respect to the 
optimal K-truncation. This is the case we mainly focus on in the sequel. The main theoretical 
results are given in subsection 2.3.1, while the main numerical results are given in section 3. ii) 
|^*'*| < Kopt(e) < \k^^^\. Then, /, is simply-resonant with respect to the optimal K-truncation. 
One such example was dealt with in the numerical study [|20|1 . Further theoretical analysis of 
this case is made in subsection 2.3.2. iii) Kopt(e) < \k^^^\ < \k'^^\. Then, /« is non-resonant with 
respect to the optimal K-truncation. Since Kop, decreases as e increases, for fixed + \k^^^\ 
this inequality always occurs if e > ei, where ei is a threshold depending on k^^\ k^^K The case 
61 > , where is the critical threshold for the onset of the Nekhoroshev regime, presents no 
practical interest. If, however, ei < ec, then, for all values of e in the interval ei < e < Cc 
the optimal normal form describes a true non-resonant dynamics. Note that in order to describe 
the dynamics close to a point of the action space corresponding to Diophantine frequencies 
a>l, it suffices to choose /» such that w, corresponds to a very high order rational approximation 
of dLt», i.e. the numbers (wi*, W2*, W3«) are high order finite digit approximants of the numbers 
(w'j^, Wj,, Wj^j). Then, |^*'^| + l^*^-*! becomes very large, and ei approaches very close to zero. In 
this case, for e sufficiently small, we expect the existence of a set of points of large measure 
within "Wi^^B, corresponding to Kolmogorov - Arnold - Moser tori in the neighborhood of the 
point It. However, these tori cannot fill an open domain. Thus, the diffusion in action space is 
topologically possible for (very weakly) chaotic orbits wandering through the set of KAM tori. 
However, in the absence of significant resonant chaotic layers (since no important resonances 
cross the question of whether or not the diffusion can be observed is of no practical 

interest, since its rate would be extremely slow to be of any relevance in applications. Thus, the 
non-resonant case is no further considered below. 

2.3.1. Double resonance 

As long as < Kopt{e), the point /» is doubly-resonant with respect to the optimal K- 
truncation. In this case, the normal form contains either terms independent of the angles, or 
trigonometric terms of the form exp(/A;*'* ■ 0^'°'"'), exp(/A;'^^' ■ 0*''"'"^) and their multiples in the 




(22) 



Kop,(e) = K'ropAe) 



(23) 
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exponents. Writing explicitly only the most important terms, the normalized Hamiltonian takes 
the form: 



3 3 , 

;=i >=i 

+ 6^/2 ^„,_„^(/'-'""))eXp(/(«ife('>+«2fc(2)).^(r„,„))^_ _ 

«i ,/i2eZ- 

The main feature of the Hamiltonian (l24l l is that, since in Z{J^'''''i"\ (/i*''"'"*)) there are coupling terms 
between more than one resonant angles, the normal form Z alone is non-integrable. In fact, Z 
can be decomposed into an integrable system of one degree of freedom and a non-integrable 
system of two degrees of freedom (see ||4|]). The decomposition is done by the linear canonical 
transformation (7f°'"', 0^'°"'^ cf>f""\ 0^'°'"') ^ (Jr, , Jr„ Jp, <Pr, , <Pr„ 0f) defined by 



(25) 



where m = {m\ , 1112, nii) has been defined in Eq.®. The Hamiltonian in the new variables reads 
(apart from a constant) 

h = Z(7«, , Jf, (pR, , 4>R,) + RopiiJR, , Jf, (pRx , 4>R2-4>f) (26) 



j(X opi) 
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4>R2 


= kf< 






(Pf'"^ + kfcpf"'' 


j(X opt) 

•'3 


- i-^l' / 

- ^3 


+ kfjR, 


+ iii^Jf, 




= mi4 


[ ('' opt) 

'1 


+ m2(fi 





where 



Z(7r, , Jr^, Jf, (f>R, , 4>R2) - ■ fn)JF 
3 



+e"^ Yj \Mij,{k'-pjR, + kfjR, + miJF)(kfjR, + kfjR, + mjJf) + ... (27) 
X gin,i>2(JRi^JR2^JF)ew('(ni4>Ri+n2(f'R2)) + ■■■ 

ni,n2&Z- 

and the remainder 7?op,(7«,, 7/7, (l)R^,(f>R^,(pF) is exponentially small in 1/e. Since (pF is ignor- 
able in Z, Jf is an integral under the flow of the normal form. On the other hand, the remaining 
degrees of freedom (/«,, 0r,) and iJR^,(f>R^) are coupled under the flow of Z due to the trigono- 
metric terms exp(/(ni0R| + na^R,))- The main characteristics of motion can be understood by the 
following remarks [j 



Since many different action symbols appeal' in tlie previous and in tlie subsequent analysis, it helps recalling that 
throughout the paper all action variables defined by a symbol starting with the letter / refer to non-scaled values, i.e. 
before the re-scaling of Eq.l llOt is implemented, while all action variables defined by a symbol starting with the letter 
/ have re-scaled values, according to Eq. UQl . Thus, in the domains considered below, all quantities of the form / - /», 
where I, is the selected central doubly-resonant point of interest, scale proportionally to e''^, while all actions denoted 
by a letter / exhibit no scaling with 6. Furthermore, all Hamiltonian-type functions denoted by ft, H'-''\ Z, or R, are 
expressed in re-scaled variables; only the original Hamiltonian (Eq.Q) is expressed in non-scaled action variables /. 
Finally, the quantities E' (Eg. 130) ) and Ez (Eq. j67t ) scale proportionally to e''^. 

12 



i) The constant- valued action Jf can be viewed as a parameter in the two degrees of freedom 
Hamihonian Z. Furthermore, except for the case of some very low resonances satisfying l^*^'^! < 
K' , all coefficients g„,,„, in ( l27l l are of order e'^^ or higher. Thus, the terms 

Zq{Jr^,Jr{, Jf) = (w, ■ m)Jf 

3 J 

+^1/2 ^ -MijAkf^jR^ + kfjR, + miJf)(k^pjR^ + kfjR^ + mjjp) + ... (28) 

'■.7=1 

define an 'integrable part' of Z, while the remaining terms depending on the resonant angles can 
be considered as a perturbation. 

The terms quadratic in , Jr^ in the rh.s. of (|28] | define the quadratic form 

^0,2 = - MiA^^JR. + kfjR.){kfjR, + kfjR^) (29) 

In Appendix A it is demonstrated that, due to the quasi-convexity condition assumed for the 
Hessian matrix M,j,, the quadratic form (|29] | is positive definite. Thus, the constant level curves 
of the quantity 

E' = (Zo - (w, ■ nijjp) (30) 



on the plane (7«[, 7r,), given by 

3 

2 



3 J 

£' = e'^' J] ^^M^Akf^jR, + kfjR, + miJpXkfjR, + kfjR, + mjjp) , (31) 
are ellipses centered at 

(/t'l) ■ M./t(2)) (m ■ M.;t(2)) - [k^^> ■ M,k^^>) [m ■ M^k^'^) 
" {km . M^kW) (/fc(2) ■ M./t(2)) - {km ■ M.fe(2))2 "^^ "-^^^ 

■ M./t(2)) (m ■ M,A:<") - ■ M,fc(i>) [m ■ M,k^^>) 

" {km ■ M,km){km ■ M,km)-{km . M.kmf 

(the role of the elliptic structures formed around double re sonances in the Nekhoroshev theorem 
is discussed extensively in f^). If the higher order terms in the action variables of the develop- 
ment of Eq.(l28]) are taken into account, the constant energy condition of Eq.(l30ll yields deformed 
ellipses on the plane (7r,,7r,). If 7r,„ or Jr^ + Jr,^, the slow frequencies = cjr^, 

4>R^ = (jjR^ are non-zero, and they are given by 

cjR, = (fe(''-M.^(i))(yR, -/«,„) + (^("-M.^(2)^(7«, J + ... (33) 

UR, - ■ M.^(2)-) ^j^^ _ j^j ^ ^^(2) . ^^^(2)-) ^j^^ _ j^ j ^ 

On the other hand, due to the definition (l25T l one has 

LOR, = ■ w(7<'''"">), LOR, = ■ w(y<'''"">), LOf^m- Loij'-'-""^) 
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which is vahd for any value of , Jr^, Jf) in the domain of convergence of the series (|28] |. It 
follows that all the resonant manifolds defined by relations of the form (nifc*''+n2^*^^)' - 
intersect any of the planes , corresponding to a fixed value of Jf. Using the notation 

A/r, = Jr, - Jr,„, aij = k^'^ ■ M,k^j\ /, = 1 , 2 

the intersection of one resonant manifold with the plane /rj) is a curve. In the linear ap- 
proximation, we have 

(niflii + n2fli2)A/R, + (niai2 + n2fl22)A7R, + . . . = 

The above equation defines a 'resonant line', which is the local linear approximation to a 'reso- 
nant curve'. All resonant lines (or curves) pass through the point (7r,„, /rjo)' which, therefore, 
belongs to the resonant junction defined by the wavevectors k^^\ k^^\ To each resonant curve we 
can associate a resonant strip in action space whose width is proportional to the separatrix width 
for that resonance. If, for a single pair of integers («i,«2), we only isolate the resonant terms 
^±n,,±n2^'^'*"''*"' in the normal form Z (Eq.(l24li). we obtain a simplified resonant normal 
form Zres(n,,m) Corresponding to the limiting case of a single resonance. In a strict sense, Z„5 
describes well the dynamics far from the resonant junction. However, it can also be used in order 
to obtain estimates of the resonance width along the whole resonant curve defined by the integer 
pair («!, n2)- To this end, the leading terms of Zr^ 5(„,_„,) are (apart from constants): 



J/2 



I 2 1 2 

-fliiAy^^ + ouAJr^AJr^ + -qiiAJr^ + 



(34) 



where the coeflicients g±„,,±„, satisfy the estimate 



-(l«,||/r<"|+|n2l|i'''l)<r 



(35) 



due to Eq.(|4]i. After still another transformation A7r, = tiiJR + n2JF, AJr^ - n2JR - rtiJf, 
(pR - n\(pR^ + n2^R2, Jf becomes a second integral of motion of Zres(nuni), which takes the form 



Zresinx^m) ~ ^ 



1/2 



c{Jf) - ^(aiiWu + 2fli2nin2 + a22nl)(JR - /r,o(-/f))^ 



(36) 



where c(Jf) and 7r,o(/f) are constants of the Hamiltonian flow of ( |36] |. Combining (l35T l and 
\, the separatrix width can be estimated as 



AJr 



32Ae-(l«ill*^'"l+l«2ll^"'l)'' 



(37) 



y flilHjj + 2ai2«i«2 + fl22n2 

Eq.(l37]i allows to estimate the width of a resonant strip in the direction normal to a resonant curve 
on the plane (7rj, Jr^). Using the relations A(A7r,) = n,A7R (for AJf = 0), this estimate takes 
the form 



AJ, 



R ,width 



32A{n\ + nl) 



1/2 



„-i(|«,||«:<i'l+|«2p<''l)o 
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(38) 




Figure 1: Schematic representation of the normal form and remainder dynamics in a domain of double resonance. Left 
panel: the resonant structure formed in the action plane of the vaiiables , ) by the overlapping of various resonant 
strips whose limits (pairs of parallel red lines) correspond to separatrix-like thin chaotic domains around each resonance. 
Two constant normal form energy ellipses E' = E\ and E' = Ei are also shown. Right: The front and back panels show 
the phase portraits corresponding to a surface of section (in one of the pairs (<l>Rf,jR^) or (<I>r2,Jr2)) under the normal 
form dynamics alone, for the energies E' = E\ (front panel) and E' = £2 (back panel). The blue curly arrows in both 
panels indicate the directions of a possible 'drift' motion (=slow change of the value of E' ) due to the influence of the 
remainder on dynamics. 
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The outcome of the analysis so far can be visuahzed with the help of Figure [T| (schematic). 
The left panel shows the structure of a doubly-resonant domain in the plane of the resonant action 
variables {Jr^,Jr^). The two bold ellipses correspond to the constant energy condition for two 
different values of E' , namely E' - E\ and E' - E2 with £1 > £2. Their common center is the 
point Jr^o) defined in Eq.(l32li. The three pairs of parallel red lines depict the borders of 
the separatrix-like thin chaotic layers of three resonances passing through the center Infinitely 
many such resonances exist, corresponding to different choices of integer vectors n = («i,«2); 
however, their width decreases as \n\ increases, according to Eq.(l38ll. We thus show schemati- 
cally only three resonances with a relatively low value of \n\, named by the letters 'A', 'B' and 
'C . The blue curly curves indicate a slow drift undergone by the chaotic orbits along the reso- 
nance layers, allowing for a transition from one resonance to another. This phenomenon, which 
will be addressed in detail below, is due to the influence of the remainder terms of the normal- 
ized Hamiltonian on dynamics. Here, however, we discuss first the (non-trivial) influence of the 
normal form terms on dynamics, by considering the Hamiltonian flow under the approximation 
H ^Z. Then, the following facts hold: 

- For any fixed value of E' , and a fixed section in the angles, the motion is confined on one ellipse. 

- For E' large enough {E' - Ei, outermost ellipse in the left panel of FigHJ, the various reso- 
nant strips intersect the ellipse E' = Ei at well distinct arcs, i.e. there is no resonance overlap. 
The right front panel in Fig{T] shows schematically the expected phase portrait, which can be 
obtained by evaluating an appropriate Poincare surface of section, e.g. in the variables (7r, , ^s,) 
or {Jr2,4>r2)- The dashed lines show the correspondence between the limits of various resonant 
domains depicted in the left and right panels. In particular, the intersection of each resonant 
strip in the left panel with the ellipse E' - Ex corresponds to the appearance of an associated 
island chain in the right panel. The size of islands is given essentially by the separatrix width 
estimate of Eq.(l38ll. Hence, the size of the islands decreases exponentially with the order of the 
resonance n - \n\\ + Inil- However, the main effect to note is that, since all resonant strips are 
well separated on the ellipse, the thin separatrix-like chaotic layers marking the borders of each 
of their respective island chains do not overlap. As a result the local chaos around one resonance 
is isolated from the local chaos around the other resonances. In fact, the normal form dynamics 
induces the presence of rotational KAM tori which, in this approximation {H ^ Z), completely 
obstruct the communication among the resonances. Note that a detailed study of the dynamics 



of the above type, induced by the doubly-resonant normal form, was recently presented in 02411 . 



- Far from the domain of resonance overlap, the size of the islands corresponding to each res- 
onance is nearly independent of the energy E' , as it depends essentially only on the size of the 
Fourier coefficient of the corresponding harmonics in the Hamiltonian. However, the separation 
of the islands is reduced as the energy decreases, since this separation is given essentially by 
the separation between the distinct arcs in FiglT]at which the various resonances intersect the 
ellipse corresponding to a fixed energy E' . As a result, below a critical energy E'^, significant 
resonance overlap takes place, leading to the communication of the chaotic layers of the various 
resonances and an overall increase of chaos. This is shown in the left panel of FiglUfor an ellipse 
E' - E2 < E'^, with the corresponding phase portrait shown in the right back panel. We note in 
particular the 'merging' of all three resonant domains one into the other, which produces a large 
connected chaotic domain surrounding all three island chains (and many other smaller chains, 
not visible in this scale). 
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The value of the critical energy marking the onset of large scale resonance overlap can be 
estimated as follows: Each resonant strip intersects one fixed energy ellipse on one arc segment. 
Also, Eq.(l38b can be replaced by the estimate 



(32A)^ 

where n - \ni\ + \n2\, k\ 2 - + |A:^^'|)/2, and M/, = {Umin + A'mflA)/2, with the constants 

Hmin,l^max defined as in Eq.©. The total length S res of all segments can be now estimated by 
summing, for all n, the estimate ( [39] l, namely 

5,„»<H4)!!|;,-**,,..<i^,-^.. (40, 

Mhki,2 ^ M/,fei,2Cr 

On the other hand, the total circumference of the ellipse for the energy E' is estimated as - 
jtR{E')^ where RiE') is the geometric mean of the ellipse's major and minor semi-axes. For 
R{E') one has the obvious estimate RiE') ~ {2E' /{e^^^Mf,)y^^, whence 

- <41, 



The critical energy E' - E[. can now be estimated as the value where S (£') ~ S res, implying that 
the associated ellipse is fully covered by segments of resonant strips. Thus 

Eq.(l42ll implies that E'^ is a (9(e'/^e"2*' -°") quantity. 

So far, we have neglected the role of the remainder in dynamics. In FiglT] the drift in action 
space caused by the remainder is shown schematically by the blue curly curves in both the left 
and right panels. Their significance is the following: The energy E = h corresponding to the total 
Hamiltonian h - Z + of Eq.(l26]l is an exactly preserved quantity. Thus, the doubly-resonant 
normal form energy E' as well as Jf cannot be preserved exactly, but they are approximate 
integrals, i.e. they undergo time variations bounded by an 6'(||^^'^°'"^||) quantity. In FigH] such 
variations will in general lead to a very slow change of the value of E' , i.e. a very slow drift 
of the chaotic orbits from one ellipse to another We seek to estimate the time required for the 
remainder to induce a transition between two ellipses with an energy difference of the same order 
as E'^, namely 

E'^-E\^ 0(e'/^e"i'^'^'^) (43) 

assuming that this effect can be described as a random walk in the value of E' (numerical ev- 
idence for this assumption will be provided in section 3). Let T be an average period of the 
oscillations of the resonant variables. By Eqs.(l34li and (l35T l. the estimate T ~ (eA)"'^^e"«^^*' -'^^^ 
holds, for a constant Ue/f ~ 1 marking the order of the most important resonances in (l34l . In 
consecutive steps, dE' can be either positive or negative, while its typical size is \dE'\ ~ \\Ropi\\- 
Then, after steps of a random walk (in the values of E'), we find an rms spread of these values 
given by 

AE ^ N"^\\Rop,\\ (44) 
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Figure 2: Same as in tlie left panel of Fig[T] but for a simple resonance. In this case, any other resonance crossing the 
main (guiding) resonance has an exponentially small width and acts as a 'driving' resonance for diffusion. 

Using ( l43l l and (l44l l. the number of steps required for the spread AE to become equal to E'^ - E[ 
(given by (l43l l is N ~ ee^"'ff'''--'^\\Rop,\\^^ ■ The diffusion coefficient can be estimated as 

i.e. the diffusion coefficient scales as ffie square of ffie size of the optimal remainder function. 
This relation is probed by detailed numerical experiments in section 3. 

2.3.2. Simple resonance 

When k^^^ < Kopt(e) < kP'\ h is simply-resonant with respect to the optimal K-truncation. 
In this case, the normal form contains terms either independent of the angles, or depending on 
them via trigonometric terms of the form exp(/«fe*'* ■ 0*^''"'"*), n e X* ■ Using the same notations as 
in ffie previous subsection, the transformed Hamiltonian reads: 

/i(7<'''""', ^('■°"'') = Z(7*''°'"', (^*''°'"') + R(J^'''''"\ 0<''""^) 

3 3 



= cj, ■ + e'/2 ^ ^ -Mijjf"'"^/!-"^ + ... (46) 
'■=1 j=i ^ 

+ e'/2 gniJ''"""^) exp(i(n/t('' ■ (^*''"'"') + . . . 

neZ' 

+ Ri/''-"'\ 0(''°'">) 

Repeating all steps as in the case of double resonance leads to the normal form 

1 9 1 7 

Zres = -ouAJr^ + anAJR.AJR, + -Uii^Jr^ + ... (47) 



(g„e'«^«. +g.„e-'"'^'>) + ... 
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The main difference with respect to the doubly-resonant normal form (l34l is that, the angle i^r, 
being ignorable, the action 7r, (or AJr^) is an integral of the flow of Zres, in addition to Jf. Thus, 
Zres defines an integrable Hamiltonian. A pair of constant values Jf = ci, AJr^ - C2 defines a 
straight line 

A/r, = -^C2 (48) 

flu 

which corresponds to the unique resonance wr, (7'''°'"') = A:*'' ■ a){J^''''"^) - 0. This will be called 
'main resonance' (- the 'guiding resonance' in [IqI]). In Figure |2] (schematic), the domain of the 
main resonance is delimited by two vertical thick red lines corresponding to the separatrix-like 
thin chaotic layers at the boundary of the resonance similarly to FiglU Using similar arguments 
as in the derivation of Eq.(l39]l, the separatrix width can be estimated as 

A7r,h.,w<. - ^^^.-^1^'""^ . (49) 

Under the normal form dynamics, motions are allowed only across the resonance, i.e. in the 
direction A7r, = const. In Fig|2]this is the horizontal direction. The thin strip delimited by two 
horizontal red lines corresponds to the resonance with resonant wavevector k^^\ which, since 
^(2) > K{e), is now of width exponentially small ((9(e'^^e"'^'**^''^^). Thus, it will be called a 
'secondary' resonance. 

In order to estimate the speed of diffusion as a function of the optimal remainder in this case, 
let us note first that the influence of the remainder on dynamics is to slowly change the value of 
the two approximate integrals Jf and AJr^, that would be exactly preserved under the normal 
form dynamics. In view of Eq.(l47li. the Hamiltonian (l46T l can be approximated by 
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-ouAJr^ + anAJR.AJR. + ^a22AJi + ... + 2/r, cos((^r,) ■ 



(50) 



h X (m- LOt)Jf + e'^^ 

where i) the (non-integer) vectors /c,, / = 1, 2, 3 come from the solution of the right Eqs.(l25t for 
the angles (p. "'" in terms of the angles 0r,, ^r,, and <pf, and ii) we approximate all the Fourier 
coefficients in the remainder series by their constant values fk* at the points A7r, = A7r, = 
(weset/R, =/i»forfe = A<i>). 

The latter approximation is sufficient for estimates regarding the speed of diffusion. The key 
remark is that for all the coefficients fk, the bound < ||/?op/|| holds, while, for the leading 
Fourier term exp(/A;/ ■ 0^''"'"*) in the remainder we have ~ ||/?op/ll- In fact, we typically find 
that the size of the leading term is larger from the size of the remaining terms by several orders 
of magnitude, since this term contains a repeated product of small divisors of the form ki ■ (see 
Appendix A). Furthermore, using an analysis as in [17], we readily find = (1 - d)K„p,, where 
< c/ < 1 is a so-called (in iITtIi ) 'delay' constant. We note in passing that the Fourier terms of 
the form exp(/A;; ■ 0^'^°'"^) are called 'resonant' in 149.1 . The value of the diffusion coefficient can 
now be estimated by applying the heuristic theory of Chirikov (|@], see also 111] and 10]) in the 
Hamiltonian model ( fSOl l. The estimate 

D ~ -^\fu,A^M\K,\f (51) 
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holds, where Qc = e''l\fl'^, T = ln(32e/w)/nc is an average period of motion within the main 
resonance separatrix-hke thin chaotic layer, of width w, A is the Melnikov function with argument 
\ki\ (see Appendix B of 112 ill ), the vector /c/ being defined by the relation a:/j(^r, +ki^2<Pr2+i(i,3'Pf - 



The estimate A(|a-/|) ~ 8;^|A:/|e "^"'^^^ holds. In view of Eq.(l25]) however, we have that 
\Ki\ = 6>((1 - d)K^°P'^ l\k^^\). Since Ko,,, ~ e"'^'* (see Appendix B), and \\R„p,\\ ~ g-'^^'""", it 
follows that Ad/f/l) ~ e"^^'*||/?opfll'', for an exponent b > Q. Putting these estimates together, we 
finally arrive at a steeper dependence of the diffusion coefficient D on the optimal remainder 
||^op,|| in the case of simple resonance than in the case of double resonance, namely: 

2Q2r 

Regarding now the precise value of b, it is hardly tractable to determine this on the basis 
exclusively of the behavior of the Melnikov integrals discussed above. We note, however, that 
the quantity A(ki) yields the size of the 'splitting'5 of the separatrix of the main (guiding) reso- 
nance due to the effects of the leading term in the remainder function. The relation between the 
separatrix splitting and the size of the optimal remainder has been examined in ifsoll and later in 
1149]. In the latter work, the estimate S ~ jj.^^^ was predicted and probed by numerical exper- 
iments, where ju (in the notation of ll49ll ) is the effective size of the perturbation to the normal 
form pendulum dynamics caused by the remainder Setting thus fi ~ \\Ropt\\ suggests the scaling 
A(ki) ~ S ~ \\Ropt\\^^^, whereby the constant b can be estimated as b ^ 1/2. Hence (in view of 
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D ~ \\Rop,f 

in simply resonant domains. 

Despite the heuristic character of the above derivation, it seems that the value b ^ 1 /2 is 
supported by the results of numerical experiments. In particular, in [20] the diffusion coefficient 
D along a simple resonance was compared directly to the size of the optimal normal form re- 
mainder. It was found that D oc \\Rop,\\^''^^ , essentially confirming that p =2(1 +b) =^ 3. We point 
out, however, that in [41J a different exponent was found p ^ 2.56 regarding the same resonance 
as in ll20ll . while it was found that p -2.1 in the case of a very low order simple resonance (with 
< K'), which is not discussed in our present work. These exponents, on the other hand, 
depend on the chosen definition of the numerical measure used to estimate both S and ||^opf||. 
Thus, a detailed quantitative comparison of the works cited above is left as on open problem for 
future study. 



3. Numerical results 

In our numerical work we employ the same Hamiltonian model of three degrees of freedom 
as in Il22[l32ll38ll33ll . The Hamiltonian reads: 



-I- /3 -I- 



4 + cos 4>\ + cos tp2 + COS 03 



(53) 



This model has a particularly simple, yet sufficient for our purpose, structure, allowing to probe 
numerically all steps of the previous section. In particular: 
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3.1. Analyticity and convexity 

The function ( l53T l is polynomial in the action variables, thus it is analytic in any complex 
extension of I - 'P? . On the other hand, the domain of analyticity in the angle variables was 
examined in ll20il . It was found that analyticity can be established in a set Re{(j)i) e T, Im((f>j) < cr, 
i = 1,2,3, for a positive constant cr estimated semi-analytically as cr ^ 0.82. Accordingly, the 
coefficients of the Fourier development 

= y y y hk expiik ■ (h) (54) 

4 + cos 01 + cos 02+ cos 03 " 

A:i=-oo A:2=-oo A:3=-oo 

where k = (ki,k2,k3), = (0i,02,03), decay exponentially. The distance of the nearest sin- 
gularity, with respect to each of the angles 0,, from the real axis is given by the solution of 
cos(p - -4/3, or = TT + 0.795365/. Thus, the following bound holds: 

\hk\ < A exp(-|A:|o-), A ^ 0.05, cr = 0.795365 . (55) 

As regards convexity, for all /, e J the matrix M» has a particularly simple structure, since 
we have Mn, - M22* - 1, and M,j, - for all other /, j. Thus there are two positive eigenvalues 
equal to unity and one equal to zero, while jL/,„,„ = i = 1 . 

The constant energy condition E - {I^ + l\)l2 + 73 defines a paraboloid in the action space. 
The resonant manifolds are planes, since wi - I\, 102 - h, <y3 = 1, whereby the resonance 
conditions 

kiLOi + k2<jL>2 + k^ciJi — kill + ^2^2 + ^3=0 (56) 

for all k = (ki,k2,kj,) define planes normal to the (hjh) plane. It follows that, when projected 
to the ijijh) plane, the intersections of all resonant manifolds with a surface of constant energy 
of the unperturbed problem yield a set of straight lines. This greatly facilitates the numerical 
study, since all diffusing orbits in the perturbed problem follow piecewise straight paths nearly 
parallel to one or more resonant lines of the unperturbed problem, while the orbits can only 
change direction by approaching close to resonance junctions. Examples of diffusion of this type 
along a simple resonance where studied in 1,38,1 , while the case of consecutive encounters with 
doubly-resonant domains was examined in a mapping model ll33l] variant of the Hamiltonian 
model ( l53T l. 

3.2. Normal form construction and optimal remainder 

The connection between the size of the optimal remainder ||/?oprll and the diffusion coefficient 
D in a case of simple resonance was the main subject of a previous study |20]. Following the 
same terminology and notations as in section 2 above, the point 7, in the normal form construction 
in i20tl was chosen as (I\*,l2*,h*) - (0.31,0.155,1). For this point we have (viz. Eq.®) 
k^^^ = (1,-2,0), = (100,0,-31), = (31,155,100). The optimal truncation order in all 
calculations of 1I20I1 varied from Kopt(e) = 18 to Kopi{e) - 39 (depending on the value of e in the 
range considered). Thus, in all cases we have < Kopt{e) < \k^^\ that is the so-chosen point 
7» was found to be simply resonant with respect to the optimal K-truncation. Following Fig. 5 of 
ll20ll it was then found by numerical fitting that the diffusion coefficient D scales with the optimal 
remainder as D oc HTJoprlP^*^. A theoretical justification for this 'steepening' of the power-law 
with respect to the exponent p =^2 holding in double resonances was given in subsection 2.3.2. 
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In order to probe now the dependence of D on ||i?op/ll in the case of a double resonance, in the 
sequel we focus our numerical study on a different point of £), namely {Iu,h*,h*) - (0.4,0.2, 1). 
The basic resonant wavevectors are 



= (1,-2,0), = (2, 1,-1), implying m = (2, 1,5). (57) 

The Hamiltonian normalization is carried out as exposed in subsection 2.2. The interval of values 
of e considered is 0.001 < e < 0.02 which, according to llssll is below the critical value for the 
onset of the 'Nekhoroshev regime' (e^. ^ 0.03). Furthermore, it will be shown below that for all 
values of e in the above interval the optimal Fourier-truncation order Kopt turns to be much larger 
than K = 4. On the other hand, for the basic wavevectors we have - 3, \k^^^\ = 4. Thus, for 
all considered values of e one has < \k^^^\ < Kopt(e), that is, the point /» is doubly resonant 
with respect to any of the optimal K-truncations considered in the sequel. 

Due to Eq. (fT2l) . the constant K' in terms of which book-keeping is implemented changes with 
e. However, one notices that, because of the logarithmic dependence of K' on e, in the largest 
part of the interval 0.001 < e < 0.02, where we focus, one has a constant value K' = 3, while one 
has K' - 2 only close to the upper limit e = 0.02 and K' = 4 close to the lower limit e = 0.001. 
For simplicity, we thus fixed the value of K' as K' - 3 in all normal form computations. Doing 
so, computer memory limitations restrict all computed expansions to a maximum order r,„a^- =17 
in the book-keeping parameter A, or maximum order \k\„,a^ - UK' - 1 = 50 in Fourier space. 
In fact, for e > 0.005 we perform at most 14 normalization steps, so that the remainder contains 
terms of at least three consecutive orders in A, namely r - 15, 16 and 17. As explained below, 
this allows us to perform some numerical tests regarding the convergence of the remainder series 
when the optimal normalization order is as high as ropt-14- (or K„pi - 42). On the other hand, 
for e < 0.005 we allow for one more normalization step (r = 15) in order to get as close as 
possible to the optimal order, which, as shown below for e < 0.004 is larger than 14. Thus, for 
the calculation of the corresponding remainder value at this order (r = 15) we necessarily have 
to rely on the sum of only two rather than three or more consecutive terms. 

Writing the truncated (at order 17) remainder function as: 

17 

7?W(/'-',0('-')<i7 = 2 R':;'\f\(P^'''>) , (58) 

s=r+l 

where R^'\ j^''\4>^''^) are the terms of order s in the book-keeping parameter A allows us to probe 
numerically the convergence of the remainder function within any chosen domain "Wi^^b in action 
space. To this end, at any normalization order r, let us consider a disk (7*'')^ + (-/2'^V ^ in the 
space of the transformed action variables (we neglect the action /3 which, in the particular case 
of the Hamiltonian ( l53l l. is dummy, i.e. it does not appear in any higher order term of either the 
normal form or the remainder). This is a deformed disk also in the old canonical variables J\,J2, 
limited by a boundary given approximately by e{J^ + J^) - {I\ - /i»)^ + (I2 - h*f - ^P^ (cf. 
the action re-scaling given by Eq.lfTOt). For the Hamiltonian ( l53T l one can readily check that all 
the terms in Ft"s\j^''\ cp^''^) are trigonometric polynomials of maximum degree K' s - 1 = 3s - 1 
whose coefficients are polynomial of maximum degree s - 1 in the actions, namely 



3.S-1 



( s-\ q \ 



\kl=0\c]=l) /=() 



exp(/fc ■ 0*''') . (59) 
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Furthermore, in the disk ^/ ,,1/2^ the obvious bound 



sup Kynvrr'l- ' (60) 



'W 



holds. We thus define the norm 



1*1=0 f/=0 /=0 \ H I 

in view of which a numerical estimate of the size of the remainder within ^i, .p can be obtained. 
In fact, by calculating the truncated sums 

.!=r+ 1 

for any fixed choice of p, where p takes all values p - r + 1, r + 2, . . . , 17, we can have a clear 
numerical indication of whether the remainder function was calculated up to a sufficiently high 
order for convergence to have been practically reached. 

The maximum value of p for which the series \\R^''\j'^''\ <;A*'^^)||<co,nv, 1/2 converges absolutely 
sets the size of the doubly -resonant domain B - e'^^Pmax (in non-scaled variables) where the 
normal form calculations are valid. In practice, we are interested in the diffusion of orbits with 
initial conditions inside this domain. In particular, in subsection 3.2 we will consider orbits 
starting on the circle po - 0.27. All our numerical orbits are studied up to a time in which their 
distance from the center of the double resonance changes significantly less than Ap = 10"' (see 
below). Variations of this order at maximum are found when we measure p either in the original 
canonical action variables or in the variables after the optimal canonical transformation. Thus, 
for all the orbits we can set a safe outer boundary p < p\, - 0.4 within which they are well 
confined. We then verify numerically that this domain belongs to the analyticity domain of the 
various transformations employed in the form of series (of the new variables in terms of the old 
variables or vice versa). This check is made by finding whether the Fourier coefficients of the 
series exhibit an exponential decay. An example is given in Fig|3] We consider the Fourier series 
yielding the new transformed canonical action as a function of the old canonical variables, 
for e - 0.01, at the normalization orders r = 4, 8 and 11. Writing this as a series 

we define the coefficients 

s,„„,-(i) .si/2 .S2/2 



h,k2, |/t,| + |lt2|=|/t|si,S2=0 



Figure |3] shows the coefficients G|A;|*'^' for e = 0.01, pi, - 0.4, and r = 4, 8 and 11. We observe 
that all three curves exhibit a tail showing exponential decay of the Fourier coefficients. However, 
it is remarkable that the asymptotic exponential slope seems to change only marginally. Instead, 

23 




* 

* 



5 10 15 20 25 
|k| 

Figure 3: The logarithm of the quantity g|^' (see text) as a function of the Fourier order \k\, for e = 0.01, at the 
normalization orders r = 4, r = 8, and r = 1 1 (lower, middle and upper set of points respectively). All three curves 
exhibit an exponential decay for large \k\, with nearly the same asymptotic law. The straight line has inclination cr = -0.8. 

the main change, as r increases, regards that formation of a 'plateau' of Fourier coefficients 
of nearly constant size formed for small \k\. Namely, the width of the plateau increases as r 
increases. It is remarkable that the asymptotic tail laws for all r appear to follow an exponential 
decay with the same constant cr ^ 0.8, i.e. with nearly the same value as the constant appearing 
in the analyticity condition of the original Hamiltonian (cf. Eq.dSSll). This effect shows that, 
while in the usual proofs of the Nekhoroshev theorem one requires a reduction of the analyticity 
domain at every normalization step, i.e. one considers bounds of the form G\k\^''^ < A*''^e"°"'''^' 
with (Tr < (Tr-i < ... < cTj, in practice the dependence of the coefficients G|^| on \k\ is more 
complicated than a simple exponential decay law. In fact, the constants o",- reflect an average 
exponential slope that compensates between the plateau, for small \k\, and the exponential tail, 
for large \k\. Namely, as the width of the plateau increases with r, one obtains smaller and smaller 
values of the average exponential decay constant cr^. 

Fig|4^ shows now an example of the behavior of the truncated remainder function for p = po 
and e = 0.01. The upper curve shows the value of \\R^''\ J^''\ 0*'^')||<p,'vi/, tp at the normalization 
order r = 6 as a function of p fox p - 7, ...17. Clearly, after p = 9 the cumulative sum ( |62] | shows 
no further substantial variation, which indicates that the remainder series converges after three 
consecutive terms p - 7,E and 9 (this is verified also by computing numerically a convergence 
criterion like d' Alembert's criterion). The lower and middle curves show now the same effect for 
the normalization orders r = 1 1 and r - 14 respectively. Note that the three consecutive trunca- 
tion orders p = 15,16 and 17 allowed for the computation of the remainder at the normalization 
order r = 14 are essentially sufficient to demonstrate the convergence of the remainder Hence, 
0*''')||<i7_^^ represents a good numerical estimator of the value of the remainder 
series for normalization orders up to r = 14. However, the main effect to note is that the estimated 
remainder value \\R^'^\J^'''\(f>^'^^)\\<\i,^^ found for r - 14 is larger than the one for r = 11, 
implying that the optimal normaUzation order ropt is below r = 14. Figll}? shows, precisely, the 
asymptotic character of the above normalization, showing ||/?*'^'(/*'^*, 0*'^')||<i7,'w, 1/2 against the 
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Figure 4: (a) The value of the remainder norm ||/?''^'||<„'u/ , „ as a function of the truncation order p when e = 0.01, 
po = 0.27, and the normalization orders are r = 6 (upper curve), r = 1 1 (lower curve) and f = 14 (middle curve), (b) The 
value of l|S*'^*ll<i7,'iV;,^u ^ function of r for different values of e. For e = 0.004 and e = 0.003, the dashed curves after 
the order r = 15 are found by quadratic extrapolation. No attempt to extrapolate was made for e = 0.002 and e = 0.001. 
(c) The optimal normalization order Top, as a function of e together with a power-law best fitting curve. 



normalization order r for various values of e as indicated in the figure. For all values down to 
e - 0.005 we now observe the asymptotic behavior, namely the size of the remainder initially 
decreases as r increases, giving the impression that the normalization might be a convergent pro- 
cedure. However, this trend is reversed after an optimal order ropi, where the remainder reaches 
its minimum value, while, for r > ropt the remainder increases with r and eventually goes to 
infinity. We also observe that for e < 0.004 the optimal order is beyond r - 15. However, for 
e - 0.004 and e - 0.003, the computed remainder values are close to the minimum. The dashed 
extensions of the numerical curves shown in Fig|4j) correspond to an extrapolation obtained by 
quadratic fitting of the available numerical points near the corresponding minima. Using this 
extrapolation, we obtain an estimate of the optimal remainder size for the values e - 0.004 and 
e - 0.003, that will be used in some calculations below. On the other hand, for e - 0.002 and 
e - 0.001, even using the extrapolation we find that the optimal normalization is beyond any 
reliable possibility to estimate given our computing limitations. 
As discussed above, the estimate ropt 

holds LMl, i.e. r^pj is expected to be a de- 
creasing function of e. FigHJ; shows the numerical estimate for rgpt as a function of e from the 
points of minima of Fig|35. The blue curve is a power-law fitting, yielding the exponent 0.52, 
i.e. very close to the one predicted by theory. 

Since the value ||.Ro;,f|| - \\R^''''''"\J^''\ (p^'''')\\<u,^v, 1/2 depends on e, from the above procedure 
we obtain numerically pairs of values (e, \\Ropi\\i^))- In subsection 3.4 below, we will numerically 
calculate the value of the diffusion coefficient D for each one of the selected values of e, thus 
allowing for a probe of the dependence of D on ||i?op(|| and a comparison with the results of 
subsection 2.3.1. 



3.3. Resonant structure 

The resonant structure in the action space (around /») can be visualized by employing the 
method of the FLI map as in 1221 13811 . We recall that the Fast Lyapunov indicator (FLI) is a 
numerical indicator of chaos, defined for one orbit by 



FLI = logio 1^(01 
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Figure 5: FLI map in the action space (surface of section (/i, /2) of tlie Hamiltonian l53l for ip^ = 0, + \ip2\ = 0, around 
the doubly-resonant point (I\,,h,) = (0.4, 0.2) for (a) e = 0.001, (b) e = 0.005, (c) e = 0.015. The color scale represents 
the computed value of the FLI (see text) in the intervals 2 < FLI < 3 (magenta, most ordered), 3 < FLI < 3.5 (blue), 
3.5 < FLI < 4 (green), 4 < FLI < 5 (orange), 5 < FLI (yellow, most chaotic). 



where f(f) is a deviation vector, i.e. in our case ^(f) = (A(f>i(t), A(f}2(t), A(f>i(t), AIi(t), Ahit), Al^it)) 
found after solving the variational equations of motion up to the time t from some initial condi- 
tions f (0). By properly choosing a threshold value FLIo ~ logjy f, orbits with FLI < FLIq are 
characterized as regular, and those with FLI > FLI() as chaotic. Furthermore, a convenient use 
of the FLI in the visualization of the Arnold web is found by producing FLI color maps ll22ll . 
Considering a grid of initial conditions in the action space, we assign to each initial condition a 
color corresponding to the FLI value found for the resulting orbit integrated up to a sufficiently 
long time (of the order 100 - 1000 periods). This allows for illustrating the resonant structure in 
action space, as shown in Fig|5] which is an FLI map in an action domain including our chosen 
doubly-resonant point (/i»,/2«) = (0.4,0.2) for three different values of e. In all three panels, 
there are resonances projecting on (/i , h) as single yellow or orange thick lines, while other res- 
onances project as strips with a green or blue interior zone delimited by pairs of nearly parallel 
yellow or red lines. As explained in |20], this difference is only due to the particular choice of 
surface of section (^3 = 0, |0i| -1- \4>2\ < 0.1, similar to |38]). Namely, the yellow lines marking 
all resonances represent the intersection of the thin separatrix-like chaotic layers formed around 
each resonance with the chosen surface of section. This produces a pair of nearly parallel yellow 
or orange lines for any resonance (of the form k ■ cj - Q) whose leading Fourier coefficient hk of 
the resonant term exp(ik ■ 0) in the original Hamiltonian expansion has a negative real part, while 
it produces a single yellow or orange thick line if Re{hk) is positive. In the latter case, the domain 
of regular orbits inside the resonance has no projection on the chosen surface of section, while 
in the former case it projects as a strip of green or blue color. 

When e - 0.001 (Fig|5t), we easily distinguish four main resonances passing through 
{Iu,h*) - (0.4,0.2). The biggest resonant domain (green, from bottom left to top right) cor- 
responds to the resonance (jJ\ - 2u)2 - 0, whose corresponding wave-vector is the basic resonant 
wavevector Similarly, the single yellow-red thick line going from bottom right to top left 
is the resonance 2wi -H W2 - =0, whose corresponding (also basic) wavevector is k^^\ We 
also clearly distinguish two resonances of order - 5, namely (jJ\ + 3t02 - W3 = (blue), 
and Swi - W2 - W3 = (green). Many other higher order resonances cross the central doubly- 
resonant point (/it,/2*) = (0.4,0.2), denoted hereafter by O, but they are not so visible in the 
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scale of FiglS^. 

The resonant strips of all previous resonances join each other forming a domain of double 
resonance around O. The extent of this domain can be determined roughly by drawing concentric 
circles around the point O. Such circles correspond to nearly constant normal form energy values, 
as can be seen by noting that, for the particular Hamiltonian function ( |53] ), the coefficients aij 
of Eq.(l34li have the values an = 5 022 = 5, and = 021 = 0. Applying Eqs. ( 12913 1I32| | 
for the particular resonant wavevectors given by (l5Ti . the doubly-resonant normal form of the 
Hamiltonian ( l53l l expressed in resonant variables takes the form 

Z^c,+6Jf + [A + (■^«2 + Jpf) + 0{e) (66) 

where the (9(e) terms are trigonometric polynomials of the resonant angles (fiR^ - (f>\ - 2(p2, 
(pR^ - 2(pi + (f>2- <f>3, while c, - e'^^0.28186... is a constant which appears only in the numerical 
values of the quantity 

Ez^Z- 6Jf (67) 

called, hereafter, the normal form energy (Ez differs from the quantity £" defined in Eq.(l30]l 
only by the constant c,). We note that the estimate 0(e) for the trigonometric terms in Z follows 
from the estimate (l35T l for the size of the corresponding Fourier coefficients, taking into account 
that ~ g-"^!*'"'! ^ according to Eq. (fT2l i. Since the angle (pf is ignorable in the 

hamiltonian ( |66] |. Jf is an integral of the flow of Z. Furthermore, since for the particular choice 
of Hamiltonian model ( |53] ) the action 73 is dummy, implying that 73 can be assigned any arbitrary 
value without affecting the dynamical evolution of any other canonical variable, we can always 
choose the value of 73 so that Jf = 0. Then, the normal form energy condition Ez - const 
implies 

^p'-5(Jl+Ji) (68) 
where p is a (9(1) quantity. Transforming to the original non-scaled action variables we also find 

(e^'^pf ^ (7, - OAf + (72 - 0.2)2 ^^g^ 

whereby e'^^p is interpreted as the radius of a circle, around O, corresponding to a constant 
normal form energy condition. It follows that the set of all possible normal form energy values 
are represented on the {I\,h) plane as a set of concentric circles around O. Three such circles 
are drawn in FigjS^, corresponding to e = 0.001 and p\ - 0.31 (outer circle), p2 = 0.27 (middle 
circle), and p3 - 0.25 (inner circle). Their main difference concerns the degree of resonance 
overlapping in each case. Namely, for a value of Ez - 0.0104, corresponding to the outer circle 
pi = 0.31, the main visible resonances of Fig.(|5t) intersect the circle in some arcs only, while 
the remaining parts of the circle lie in the regular (non-resonant) domain. In the latter parts, 
the normal form dynamics alone would imply the existence of a set of Kolmogorov-Arnold- 
Moser invariant tori of large measure. On the contrary, in the inner circle, corresponding to 
Ez = 0.0099, all resonances essentially overlap, producing a strongly chaotic domain. The 
middle circle corresponds to Ez - 0.01007, which is close to the critical energy below which 
resonance overlapping dominates the dynamics. 

The remaining panels of Fig|5]show what happens when e is increased by a factor 5 (e = 
0.005, Fig. I5J)), or 15 (e = 0.015, Fig. |5};) with respect to Fig|5^. A main feature to notice is 
that, by increasing e, many more resonances 'show up' in the FLI map. Furthermore, the size 

27 



of all resonant domains grows proportionally to e'^^, as verified in Fig|5] where by augmenting 
the scale in panels (b) and (c) by a factor Vs and VTs respectively with respect to panel (a), the 
widths of all resonant strips passing through O remain essentially unaltered in all three panels. 
Thus, the only essential change is the increase of chaos as e increases. Namely, we see that the 
chaotic layers delimiting the borders of each resonance become thicker as e increases. This also 
increases the resonance overlapping locally, close to the points of resonance crossings. 

Focusing, now, on one value e - 0.008, Figure |6] shows in detail the implications of normal 
form dynamics in the two regimes when there is no resonance overlap {Ez - 0.0306, Figs. 
|6^,c), or when there is substantial resonance overlap {Ez - 0.029, Figs|6]3,d). The upper panels 
correspond to FLI maps as in Fig.©. Here, however, instead of the action variables (/i,/2) 
we use the resonant re-scaled actions defined as in Eq. (l25b . where, for each point 

in the action space of the original variables we compute the values of the transformed actions 
J. 1,2,3 by the composition of the Lie canonical transformations resulting from the 

computer- algebraic program calculating the optimal normal form Z. Since the same program 
renders also the algebraic form of Z, we use this expression to derive the Hamiltonian equations 
of motion of the normal form alone, namely - dZjdJR^, <pR^ = dZ/dJu^, Jr^ - -dZjd^R^, 
Jr^ - dZld(pR^, (pf - dZjdJf, while we set Jf - const - 0. For each value of Ez, we then 
compute numerical orbits under the normal form dynamics alone via the previous equations. 
Finally we plot a convenient surface of section of the normal form flow, taken by the condition 
mod{(pR^ - 2(pR^,2n) - mod{5<p2 - <p3,2jT) = 0. These sections are shown in Figs|6};,d, for the 
normal form energy values Ez = 0.0306 and Ez = 0.0290 respectively. The corresponding 
circles, through Eq.(l68t. are shown in panels (a) and (b), superposed to the color background 
yielding the FLI map in the resonant action variables for e = 0.008. The main feature of this 
plot is the exact correspondence between the values of 7r, where each resonance intersects the 
circle corresponding to Ez = const in panels (a) and (c), and the projection of these values to 
thin chaotic layers delimiting the same resonance in the corresponding surface of section. In 
fact, inside each resonance we have regular orbits corresponding to islands of stability on the 
surface of section. Furthermore, while at the normal form energy value Ez = 0.0306 there are 
many rotational KAM tori separating these resonances, at the value Ez - 0.029 these tori are 
destroyed and substantial resonance overlap takes place. This fact leads to the creation of a 
connected chaotic domain surrounding all main resonances in the surface of section of Fig|6}l. 
This, in turn, implies that under the normal form dynamics alone no communication is allowed 
from one resonance to the other for the normal form energy value Ez = 0.0306 (which in this 
approximation remains constant in time), while such communication is possible throughout the 
whole connected chaotic domain for Ez = 0.029. In fact, the phase portrait of Fig|6}l renders 
visually clear that chaos is rather strong in this case. However, as emphasized in section 2, this 
fact has no consequences regarding the possibility of long excursions in the action space, since 
all motions in this approximation would be bounded on circles like those of Figs|6t,b. On the 
contrary, such excursions are only possible due to the eff'ect of the remainder, which causes the 
chaotic orbits to slowly 'drift' from circle to circle as the value of Ez changes slowly in time. To 
this we now turn our attention. 

3.4. Visualization of Arnold diffusion in doubly-resonant normal form variables 

The main eff'ect, of local diffusion within the doubly-resonant domain, can now be demon- 
strated with the help of Figure |7] The time evolution of one chaotic orbit is shown in this figure, 
as the orbit moves within the doubly-resonant domain along some of the main intersecting reso- 
nances. In this example as well we take e - 0.008 (as in Fig|6]l. 
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Figure 6: (a) and (b) The FLI map in the plane (Jl^^ , Jr^ ) defined as in Egs. )251 for the Hamiltonian )53K for the same 
surface of section and using the same color scale as in Fig.^5). The circle in (a) corresponds to the constant normal 
form energy value Ez = 0.0306 in (a) and Ez = 0.029 in (b). The phase portraits of the normal form dynamics 
for the values (c) Ez = 0.0306 and (d) Ez = 0.029. The plotted surfaces of sections are (<f>l^^ , Jr^ ) whenever the 
quantity (j>R^ - 2^R| = 502 ~ 03 (where all symbols denote the new canonical coordinates and momenta after the optimal 
Lie normalization) crosses a multiple value of 2n. The main resonances are identified as: li>\ - 2oj2 = (vertical), 
2wi + W2 - W3 =0 (horizontal), oj[ + 3co2 - 0J3 = (bottom left to top light diagonal ), 3oji - 012 - coj =0 (top left to 
bottom right diagonal ). 
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Figure 7: Visualization of Arnold dilfusion in appropriate variables of the doubly-resonant normal form, for a numerical 
orbit in the Hamiltonian (53) for e = 0.008. After computing the optimal normal form, we find, via the Lie canonical 
transformations, the values of all transformed variables (?), (/), 7;r(r) and (;), (^r, (;), (^;r(r) corresponding to 
particular values of the old variables Ji(t), J2(t), J3(t) and 0i(O, i^2(0. 03(0 stored at many different times t within an 
interval < f < 1.5 X lO' along the numerical ran. Using the numerical values of the computed transformed variables, 
(a) shows the variation of the normal form energy Ez{t) as a function of / in the intervals < ? < 3 X 10** (blue), 
3 X 10* < ? < lO' (red), and lO' < f < 1.5 X lO' (green). The initial and final values are equal to Ez{t = 0) = Ez(t = 
1.5 X 10') = 0.0306, while the minimum value, occurring around ? = 8 X 10** is Ez = 0.029. (b) The evolution of the 
orbit in the action space (7r| , ), using the same colors as in (a) for the corresponding time intervals. In the first time 
interval (blue), the orbit wanders in the thin chaotic layer of the resonance a>i + 'iaii - a/} = 0. In the second time interval 
(red) it jumps first to the resonance 3(ij[ - - = 0, and then to the resonance toi - 2(t)2 = 0. In the third time interval 
(green) the orbit recedes from the doubly-resonant domain along the resonance toi - 2tL)2 = 0. (c) 3D plot in the variables 
(<pRj , Jri , Ez), visualizing Arnold diffusion for the same orbit. Taking 20 equidistant values of Ezj, i = 1, 2, ... 20 in the 
interval 0.029 < Ez < 0.0306, we first find the times in the interval < f < 9 X 10** when the normal form energy 
value Ez(t) of the numerical orbit approaches closest to the values Ezj. For each i, starting with the momentary values 
of all resonant variables at , we then compute 1000 Poincare consequents of the normal form flow on the same section 
as in Figs|6j;,d. The same procedure is repeated in a second interval 9 X 10** < f < 1.5 X lO'. As a net result, the orbit at 
the beginning and end of the calculation is found on the same section (corresponding to Ez = 0.0306), but in a different 
resonant layer, having by-passed the barriers (invariant tori of the normal form dynamics) via a third dimension (here 
parameterized by the time- varying value of Ez)- 
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In Fig|2l the evolution of the orbit is shown for a total time t - 1 .5 x 10^. The optimal normal 
form for e = 0.008 has also been computed, whose optimal normalization order is r^p, -12, 
corresponding to an optimal Fourier order K„p, - 36. Since the corresponding Lie generating 
functions are known, we compute, via the composition of Lie canonical transformations, the 
values of all transformed variables 7^,(0, 7R,(f), 7f(f) and (pn^it), (f>R^(t), (f>f{t) corresponding to 
particular values of the old variables Ji(t), J2(t), Ji{t) and (p\{t), (p2(t), <p3(t) stored at many differ- 
ent times during the numerical run, i.e. as t varies within the interval < f < 1.5 x 10''. Finally, 
since the exact algebraic expression for the normal form Z is known, we compute the precise 
numerical value of the normal form energy Ezit) at the same times. 

FigEh shows the variation of the normal form energy Ez{t) as a function of the time t in the 
intervals < f < 3 x 10** (blue), 3 x 10'^ < r < 10'' (red), and 10'' < f < 1.5 x 10'' (green). 
The final time is such that the initial and final values of Ez are equal, namely Ez{t = 0) = 
Ez{t = 1.5 X 10^) = 0.0306. On the other hand, as Ez slowly changes during the run, it acquires 
a minimum value around f = 8 x 10*^, which is Ez^min = 0.029. Such evolution corresponds 
to the process described schematically in Fig{T] (section 2). Namely, from the previous figure 
(Fig|6ll we conclude that the two extreme values of Ez acquired during the numerical run are 
such that Ez{t = 0) > Ezc while Ez^min < Ezc, where Ezc is the critical energy corresponding to a 
large scale overlapping of resonances (subsection 2.3.1). Furthermore, as we will see in the next 
subsection, the chaotic excursions of the orbits, and, consequently, time evolution of Ez, can be 
approximated by a normal diffusion process. Furthermore, the fastest evolution takes place in the 
intervals < r < 10^ and 1.3 x 10^ < f < 1.4 x 10^, in both of which the total variation of Ez is 
of the order of 10"^, or a 'per step' variation of the order of AEz ~ 10 It should be stressed 
that these extremely small variations are possible to unravel numerically only because we use the 
new canonical variables deduced by the normalizing sequence of Lie canonical transformations. 
When the old variables are used, instead, we find that the there are large variations (of order 
e'^^) of all quantities depending on the actions. These variations are, in fact, dominated by the 
so-called (in the Nekhoroshev theory) 'deformation' effects (which are also of order e'^^), hence 
completely covering the drift eff'ects which are much smaller in size. This feature of the optimal 
canonical transformations will be exploited in the measurement of the diffusion coefficient D as 
described in the next subsection. 

FigEJ) shows the diffusion of the orbit in the action space (/s, , /«,), using the same colors as 
in Fig|7^ for the corresponding time intervals (the background produced by the FLl map is shown 
here in gray scale). In the first time interval (blue), the orbit wanders chaotically within the thin 
chaotic layer of the resonance a>i + 3a)2 - W3 = 0. It should be stressed that this wandering has 
a random walk character, i.e. the orbit makes several reversals of its drift direction, sometimes 
approaching and other times receding from the center of the double-resonance. On average, 
however, the drift is in the inward direction (this is a statistical effect; for other initial conditions 
the average drift turns to be outwards). In the second time interval (red), the orbit jumps first to 
the domain of the resonance 3wi - W2 - ^^3 = 0. Now, however, the chaotic motion takes place 
with a relatively high speed (of order e'^^) in the direction across resonances. As a result, the orbit 
fills nearly ergodically the whole connected chaotic domain surrounding the main overlapping 
resonances, while, at the end of this time interval, the orbit is closer to the resonance wi -2w2 = 0. 
Finally, in the third time interval (green) the orbit recedes from the doubly-resonant domain (this 
is also a statistical effect) being trapped along the domain of the resonance a>i - 2a>2 = 0. In this 
way, at the time f = 1.5 x 10^, the orbit is found at about the same distance from the center as 
initially (at t = 0), but on a different resonance. 

FigEt, now, shows a 3D plot in the variables ((pu^ , Jr^ , Ez), visualizing the 'third dimension' 
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along which the Arnold diffusion progresses for the same orbit. From this plot we can clearly 
see the effect of the remainder, which can be considered as a very slow modification of the 
normal form dynamics acting on a timescale of the order of 10^ periods. The normal form 
dynamics, on the other hand, describes well the motion over shorter timescales, of the order of 
lO'^-lO^ periods. In order to show the dynamical effects happening on both timescales, we adopt 
the following numerical procedure: Taking 20 equidistant values of Ezj, i - 1,2, .. .20 in the 
interval 0.029 < Ez < 0.0306, we first find the times f, within the interval < f < 9 x 10^ (where 
the motion is, in general, in the inward direction) when the normal form energy value Ez{t) of 
the numerical orbit approaches the closest possible to the values Ezj- Then, for each /, we set 
the momentary values of all canonical variables of the numerical orbit at the time f,- as initial 
conditions via which we compute the corresponding values of all the new resonant canonical 
variables following the composition of the corresponding Lie canonical transformations. With 
these values as initial conditions, we compute 1000 Poincare consequents of the normal form 
flow alone on the same surface of section as defined in Figs|6};,d. The same procedure is repeated 
in the second interval 9 x 10*^ < f < 1.5 x 10"^, where the motion is in general in the outward 
direction. The whole set of Poincare consequents (points {(I>r,,Jr,) gathered in this way are 
plotted in the 2D sections of the parallelepiped of Fig]?};, along with the variations of the value of 
the normal form energy Ez(t) (sampled more frequently) which are shown in the third dimension. 

The details of the filling process of the various resonant chaotic layers located in the doubly- 
resonant domain are now clearly seen. In particular, we note that the chaotic orbit fills the 
whole separatrix layer of the initial resonance coi + 3a>2 -0)3 = in a timescale much shorter 
than the one required for substantial drift in the Ez direction. After a transient 'back and forth' 
motion around E^ = 0.03, the orbit then moves slowly towards the value Ez = 0.029, where all 
important resonances overlap. In the intermediate time interval (red), we clearly see the filling 
of the stochastic layers of both resonances 3dLti - cj2 - (^3 =0 and wi - 2w2 = 0, while global 
transport is allowed by the normal form dynamics from one resonance to the other. As, however, 
the remainder effect causes a new motion of the orbit outwards (i.e. towards higher values of Ez 
(green)), the orbit is eventually captured at the resonance cji - 2(l>2 - 0, and stays there until the 
end of the simulation at f = 1.5 x 10^. 

It should be emphasized that the fact that the orbit moves in the outward direction at f = 
1.5 X 10'' does not guarantee that there will be no further return inwards. In fact, we find that 
most orbits undergo several 'in-out' cycles like the one described in Fig|7] before eventually 
abandoning the doubly-resonant domain. As an estimate, for e = 0.008 we find that the number 
of cycles before a final exit from the doubly -resonant domain is of the order of 10, while the total 
time required for this effect is of the order of 10'" to 10" periods. Furthermore, the probability 
of exit along one particular resonance decreases as the order of the resonance increases. This 
is expected, since the width of resonances scales with their order |^| as ~ e^'^l^^l/^^ while the fast 
filling of the innermost chaotic domains where all the resonances overlap is nearly ergodic. 

Finally, we point out that a visualization of the diffusion process like in Fig|7}; clearly sug- 
gests that the diffusion is driven by the intersections of the asymptotic manifolds of lower- 
dimensional objects (like hyperbolic 2D tori) all along the path in which the diffusion takes 
place. However, locating such tori, and studying their manifolds is a task that cannot be accom- 
plished by the use of the Birkhoff normal form as above, on the other hand, the latter provides 
good initial conditions for a numerical search of such tori. This subject is proposed for future 
study. 
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Figure 8: Three orbits with initial conditions on the circle R(0) = e''"Po with po = 0.27, for (a) e = 0.004, (b)) e = 0.007 
and (c) e = 0.01. The black points show the orbits' consequents on the surface of section up to a time t = 10^. All three 
orbits are diffusing outwards. The circle with radius R(0) is shown in pink. 

3.5. Dependence of the diffusion coefficient on the optimal remainder 

Our final goal is to obtain numerical estimates of the value of the diffusion coefficient D as 
well as its relation to the size \\Ropi\\ of the optimal normal form remainder as e is varied in the 
interval 0.003 < e < 0.020. 

To this end we implement the following numerical procedure: For any fixed value of e, using 
the information from the FLI maps, we first select 100 initial conditions corresponding to on a 
circle defined as in Eq.d69]l, where the radius is chosen equal to 

p = Po = 0.27 . 

For such a choice of p, the corresponding circle lies inside the resonance overlap domain, en- 
suring that the short time dynamics is dominated by the doubly-resonant normal form. How- 
ever, in longer times all these orbits exhibit weakly chaotic diffusion. The complete set of ini- 
tial conditions for one orbit on the circle p - po - 0.27 are found by solving simultaneously 
for 1 1 and I2 the equation of the circle (Eq.(l68]l) as well as an equation for the initial angle 
(po = arctan[{l2- 0.2) /{I I -0.4)], where, for each initial condition, 0o is chosen by visual inspec- 
tion so as to correspond to an initial condition in the domain of each one of the main overlapping 
resonances. 

We then follow numerically these orbits for a time long enough so that the mean change of 
their radial distance from the center is large enough to allow for a reliable computation of the 
diffusion coefficient. Let R{t) = e^^^p{t) be the instantaneous value of the distance from the 
center for any such orbit. The quantity [Rit) - RiO)]^ changes as an orbit slowly drifts from one 
circle to another Figure |8] shows this effect for three orbits corresponding to the same initial 
angle 0o but for three different values of e, namely e = 0.004 (Fig[8h), £ = 0.007 (Fig[8}3), and 
e = 0.01 (Fig[8]:). The orbits are shown by the black points on the section |0i| + \(p2\ < 0.1, 
03 = 0, superposed as usually to the colored background of the FLI map. The pink circles in 
each panel are the circles R{Q) - e^^^po, where the orbits' initial conditions lie. 

Apart from an overall change of the size of the circle of initial conditions with e, a simple 
visual comparison of the three panels suffices to conclude that they imply quite different diffusion 
rates of their depicted orbits. In all three panels, the orbits (black points) are shown up to a time 
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Figure 9: The time evolution of the quantities Ez and ]f (see text), using the new transformed canonical variables (first 
and second panel), and p{() and ]f using the original canonical variables (third and fourth panel), for two chaotic orbits 
of our chosen ensemble for e = 0.008 (upper and lower row). 

f = 10**, which is quite long compared to the time needed to fill the chaotic domain along the 
circle p - po- However, when e - 0.004 (Fig[8^), the orbit's plot shows that the orbit exhibits no 
discernible transverse motion with respect to this circle, despite the fact that the orbit lies entirely 
within a rather strong chaotic domain (yellow in the FLI scale). On the other hand, when e is 
raised to e = 0.007 (FigHJ)), the orbit is observed to create a small ring around its initial circle, 
implying that the diffusion is visible in this timescale. Increasing e still further (e = 0.01, Fig|8};), 
causes now a rather fast diffusion, which leads to the orbit following clearly a preferential 'exit 
resonance', where the diffusion continues essentially as in the simple resonance case (subsection 
2.3.2). 

A key remark, now, is the following: similarly to the case of the orbit of Fig|7] whose dynam- 
ical features were possible to unravel using the new, i.e., transformed canonical variables after an 
optimal normalizing transformation, exploiting the same variables, instead of the original ones, 
allows to observe the random walk-like drift of one orbit in the action space in a much shorter 
integration time than by the use of the original variables. An example is given in Fig|9l for 
e = 0.008. We compute, via the optimal normalizing canonical transformation, a time sequence 
of the values of all the transformed canonical variables (y^''°'"'(f), 0*''°'"*(f)) from the available se- 
quences of values of the original variables J{t), (p(t) along the numerical orbits. The four panels in 
each row show the time evolution, for one chaotic orbit on the circle po = 0.27, of the quantities 
i) Ez computed in the transformed canonical variables, ii) Jf - (27i + J2 + 5Ji)/3Q computed 
in the transformed variables, iii) p(f) computed in the original canonical variables, and iv) Jf 
computed in the original variables. We note immediately the gain by passing the data through 
the optimal normalizing transformation, namely the fact that this transformation absorbs all 'de- 
formation' effects, allowing to see the very slow drift due to the weakly chaotic diffusion in a 
timescale t ~ 10^. In fact, the quantity Ez can only be computed in the transformed canonical 
variables, in which, for both orbits, it undergoes variations of the order 10""*. In comparison, the 
analog of Ez in the original variables, i.e., p(f), undergoes variations in the second digit, and the 
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Figure 10: The time evolution of the quantity a^^ (see text), in our ensemble of numerical data for e = 0.008, when 
computed by use of the new transformed canonical variables (left), or the original canonical variables (right). The 
evolution is shown up to the time 6 X 10^. 

corresponding time evolution is dominated by (9(e'^^) oscillations, which completely hide the 
slow drift process in the radial direction with respect to the central doubly resonant point. The 
comparison is even more straightforward in the variables Jf computed by the transformed and 
by the original action variables. In the former, we can clearly see the drift phenomenon for both 
orbits, which results in a slow change of the value of Jf (which is an approximate integral) at 
the fifth digit. In contrast, this phenomenon is completely hidden when Jf is computed in the 
original variables, since the corresponding plot is dominated by oscillations of at least one order 
of magnitude larger amplitude than the drift effect. 

In order, now, to measure the value of the diffusion coefficient, using the data from all 100 
orbits, we define the mean square deviation: 

J 100 

^M-—J]{yi{t)-y{t)f (70) 

i=l 

where y{t) = Y(t) - Y(Q), and Y(t) stands for any of the four quantities shown in Figj9] Plotting 
(Ty against the time t allows to estimate the diffusion coefficient. Figure [10] shows an example 
of this calculation, setting Y equal to Jf in the transformed variables (left panel), or the original 
variables (right panel). We note again that it becomes possible to observe the diffusion in a 
timescale f ~ 10^ using the ensemble of data in the transformed variables, while this time is 
quite short to reveal any linear trend of crj^ with the time f in the original variables. In fact, in 
the original variables it was possible to measure reliably the diffusion coefficient only after an 
integration time t - 10^. Furthermore, this time increases even more for smaller values of e. 

Figure [TT] shows the final result. Computing, as indicated above, the diffusion coefficients 
Df^ and Dj^ in the transformed canonical variables, for eleven different values of e as noted in 
the caption, we also use the data from FigH] whereby we obtain the optimal remainder value 
ll^opfll for the same values of e (from the minima of the curves of Figj4). We then plot Df^ and 
Dj^ against \\Ropt\\ in a log-log scale. Despite some scatter, the correlation of both independent 
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Figure 1 1 : Log-log plot of the dependence of the two estimates of the diffusion coefficient De^ (upper set of points) 
and Djp (lower set of points) on the optimal normal form remainder using numerical data from the integration 

of orbits (see text). The points con-espond to the values of e (from left to right) 0.003, 0.004, 0.005, 0.007,0.008, 0.01, 
0.012, 0.013, 0.015, 0.018, 0.020. The straight lines represent the power-law fits Xogi^iDE^) = -5 + 2.3 logio(||Rop,||) 
(upper) and logio(D£,^) = -7.1 + 2.2 log 10(||Rop,||) (lower). 

estimates of the diffusion coefficient with can be described as a power-law. The power- 

law exponents found by best-fitting are = 2.3 for the data of De^ and p - 2.2 for the data of 
Djj,. In these best fittings we excluded the two points for e = 0.003 and e - 0.004, since the value 
of the optimal remainder found by extrapolation is uncertain for these values of e. However, we 
note that the corresponding points in FiglTT]are still very close to the fitting law found by the 
remaining data. 

The exponents found in Fig|4]are not far from the theoretical estimate p - 2 derived in sec- 
tion 2 (Eq.(l45]l). However, we have made various trials to determine p via alternative definitions 
of the diffusion coefficient, and we always find estimations of p somewhat larger than 2. We ffius 
conjecture ffiat this difference from p = 2 is a real effect (not due to numerical uncertainties), 
which, however, requires a more detailed theory to interpret. On the other hand, the correspond- 
ing analysis for simple resonances (subsection 2.3.2) as well as the numerical results of [i20ll 
indicate that the steepening of the power law in simple resonances of order not smaller than K' 
is quite substantial, leading closer to p ^ 3. In the latter case, another independent example fl2| 
yields p ^ 2.5. The issue of how exactly to quantify ffie steepening of the power-law remains 
open. 

4. Conclusions 

We examined in detail the phenomenon of weak chaotic diffusion in doubly or simply res- 
onant domains of Hamiltonian systems of three degrees of freedom satisfying the necessary 
conditions for the holding of the Nekhoroshev theorem. The aim was to determine a quantita- 
tive relation between the diffusion coefficient D and the size of the optimal remainder ||/?op/|| of 
a resonant normal form constructed according to the requirements of the analytical part of the 
Nekhoroshev theorem. Our main results are ffie following; 
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1) We propose an efficient algorithm for Hamiltonian normalization, which is implemented 
as a computer algebraic program performing expansions up to a high order We explain the 
practical aspects of this algorithm, and show how it can be used in order to compute i) the 
optimal normalization order ropt as a function of the small parameter e, and ii) an estimate of 
the size of the remainder \\Ropt\\ at the order r^p,. The dependence of r^pt on e is found to be an 
inverse power-law with an exponent in agreement with theory. 

2) We construct estimates on the speed of diffusion in doubly resonant domains. To this 
end, we examine first the dynamics under the Hamiltonian flow induced by the normal form 
alone (i.e. neglecting the remainder). The role of the convexity conditions assumed for the 
original Hamiltonian is analyzed in the context of the normal form dynamics. We then discuss 
the influence of the remainder on dynamics. Estimates on the value of the diffusion coefficient D 
are quantified by considering a 'random walk' model for the slow drift of the value of the normal 
form energy due to the remainder The final prediction is a power-law estimate D ~ \\R„pi\\'' with 
p ^ 2 in doubly resonant domains. 

3) We perform detailed numerical experiments aiming to test the above predictions, em- 



ploying the same Hamiltonian model as in Il22ll as well as the 'FLI map' method. Using the 



information from the computed normalizing canonical transformations, we propose a convenient 
set of variables in which the Arnold diffusion in the doubly resonant domains is clearly visu- 
alized. Furthermore, using ensembles of chaotic orbits, we make two independent numerical 
calculations of the diffusion coefficient D for various values of e. The relation between D and 
||7?op,|| found by the two calculations is D ~ ||/?o/«lP'^ and D ~ \\Ropi\\'^'^ respectively. 

4) Finally, we make some theoretical estimates on the relation between D and in sim- 

ply resonant domains. In this case, we combine the basic theory developed in L2J together with 
estimates given in ll49ll regarding the dependence of the size of the separatrix splitting on the 
optimal normal form remainder in simply resonant domains. We are thus led to the prediction 
\\Roptf^^^''\ where b^l/2,0Tp^ 2(1 + b) ^ 3, holding for all simple resonances of order 
higher than K', where K' is defined in Eq. (fT2l i. The latter result interprets the results obtained in 
an earlier study li20ll by purely numerical means. 
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Appendix A. Quasi-convexity and normal form energy constraints 

The quadratic form ^0,2 given by Eq.(|29l) can be written as: 

^0,2 = Ur.Jr,) ■ fe*'-'^ ■ M, ■ (^(i-^*)^ • (A.l) 
38 



where is a 2x3 matrix whose first and second line are given by ik^^\ k^2^ , A:^'') and (^j^', ^2^', ^3^') 
respectively. Since the matrix M, is real symmetric, it can be writen in the form M, = X • //* • X'- , 
where - diag(fi],fj.2,fJ^3), with yU, = the eigenvalues of M», while X is an orthogonal matrix 
with columns equal to the normalized eigenvectors of M,. Using the above expression for M», 
Eq. dA.lb resumes the form 

where Y = A:*''^^ ■ X is a 2 x 3 matrix. Writing ^0,2 as ^0,2 = QJrj + ^-^^i Jri + ^-^r,' and denoting 
by yij the elements of Y, the discriminant A = 4-QP - is given by: 

A = -[(yil}'22 -}'l2}'2l)Vl/'2 + (3'11}'23 -}'l3}'2l)Vl/^3 + (}'12}'23 - }'l3}'22)V2/^3] (A.2) 

Since we have assumed (subsection 2.1) that either all three eigenvalues /i, have the same sign, 
or two of them have the same sign and one is zero, by Eq. (IA.2) we have that A < 0. That is, the 
quadratic form ^0,2 is positive definite. 
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